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4 DIMENSIONAL MAGNETIC RESONANCE IMAGING 

This application claims priority to U.S. provisional application serial No. 
60/226,667, filed on August 22, 2000, the complete disclosure of which is 
5 incorporated herein by reference. 

1. Field of the Invention : 

10 This invention relates to method and apparatus for imaging an object, such as a 

body. More particularly, the invention provides high resolution, real-time images in 4 
dimensions with little or no deterioration from motion artifact. 



15 2. Background of the Invention : 

Nuclear Magnetic Imaging (NMR) which is commonly called magnetic 
resonance imaging (MRI) entails 1 .) magnetizing a volume with a constant primary 
magnetic field in a z-direction, 2.) providing a gradient along the axis of the z-directed 

20 field to select a slice in the xy-plane, the plane perpendicular to the direction of the 
primary field, 3.) providing electromagnetic radiation resonant with the Larmor 
frequency of protons in the slice, 4.) providing a pulse of resonant electromagnetic 
radiation to flip the magnetization vector into the transverse plane or plane of the slice, 
and 5.) applying a magnetic field gradient along an axis in the xy-plane of the z-directed 

25 field with excitation at the Larmor frequency to provide phase dispersion of the NMR 
signal along the axis to encode spatial information, and 6.) recording the free induction 
decay (FID) radio emission signals following excitation, 7.) recording a plurality of such 
FIDs, each recorded following an excitation with a rotated direction of the gradient in the 
xy-plane, and 8.) reconstructing the image from the plurality of the FIDs. An integer n 

30 of FIDs each having a phase gradient that corresponds to the magnetic field gradient that 
was rotated to n unique directions in the xy-plane comprise a set along two orthogonal 
axes in phase or k-space. A two dimensional Fourier transform of the data set is used to 
reconstruct an n by n pixel image. 

MRI is of primary utility in assessing brain anatomy and pathology. But long 

35 NMR relaxation times, a parameter based on how rapidly excited nuclei relax, have 

prevented NMR from being of utility as a high resolution body imager. The most severe 
limitation of NMR technology is that for spin echo imaging n, the number of free 
induction decays ("FIDs"), a nuclear radio frequency energy emitting process, must 
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equal the number of lines in the image. A single FID occurs over approximately 0.1 
seconds. Not considering the spin/lattice relaxation time, the time for the nuclei to 
reestablish equilibrium following an RF pulse, which may be seconds, requires an 
irreducible imaging time of n times 0.1 seconds, which for 512 x 512 resolution requires 
5 approximately one minute per each two dimensional slice. This represents a multiple of 
1500 times longer that the time that would freeze organ movements and avoid image 
deterioration by motion artifact. For example, to avoid deterioration of cardiac images, 
the imaging time must not exceed 30 msec. A method for speeding NMR imaging flips 
the magnetization vector of the nuclei by less than 90 degrees onto the xy-plane, and 

10 records less FIDs. Such a method, known as the flash method, can obtain a 128 x 128 
resolution in approximately 40 seconds. Another technique used to decrease imaging 
time is to use a field gradient and dynamic phase dispersion, corresponding to rotation of 
the field gradient, during a single FID to produce imaging times typically of 50 msec. 
Both methods produce a decreased signal-to-noise ratio ("SNR") relative to spin echo 

15 methods. The magnitude of the magnetization vector which links the coil is less for the 
flash case because the vector is flipped only a few degrees into the xy-plane. The echo- 
planar technique requires shorter recording times with a concomitant increase in 
bandwidth and noise. Both methods compensate for decreased SNR by increasing the 
voxel size with a concomitant decrease in image quality. Physical limitations of these 

20 techniques render obtaining high resolution, high contrast vascular images impractical. 

SUMMARY OF THE INVENTION 



It is an object of the invention to provide high resolution multi-dimensional 
25 images of an object, such as a body, tissue, or working cardiopulmonary system. 

It is a further object of the invention to rapidly acquire the data to provide 
magnetic resonance images of a body with reduced motion artifacts. 

These and other objects of the invention are attained by providing an apparatus 
for obtaining a magnetic resonance image of a body using data acquired over the three 
30 spatial dimensions plus time, rather than acquiring data only in time at one receiving 
antenna. A preferred embodiment of the apparatus of the invention includes a radiation 
source for applying a first radiation field having a magnetic component to the body, to 
magnetize the body. The apparatus further includes a source for applying a second 
radiation field to the body, to elicit a radiation field from the body. A detector senses 
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this radiation field, and produces a signal that a reconstruction processor employs to 
create the magnetic resonance image of the body. 

A NMR image is obtained of a magnetized body from a three-dimensional map 
of the intensity variation of the NMR signal produced by each voxel of the magnetized 
5 body, and detected over a three-dimensional volume of space external to the body, herein 
referred to as the "sample space." The data is acquired over three spatial directions plus 
time. In an embodiment, the NMR signals are detected over a three dimensional detector 
array as a function of time. The NMR signals may be sampled at least at the Nyquist 
rate, i.e., at a rate that is twice the highest temporally frequency of the NMR signal and 

10 twice the highest spatial frequency the Fourier transform of the NMR image of the 
phantom. Sampling at the Nyquist rate or higher allows the spatial variations of the 
external NMR signal to be acquired. 

In an embodiment, the NMR signal at each detector as a function of time is 
processed by a method such as a Fourier transform operation to give a plurality of 

15 Fourier components each having the same frequency, an intensity and a phase angle. 
The NMR signal of each voxel at any given detector gives rise to a Fourier component 
with a unique phase angle relative to the Fourier component of any other voxel of the 
phantom at that detector. The set of Fourier components that correspond to the NMR 
signal of a given voxel over the detectors is determined. This may be achieved by using 

20 a first component having a phase angle and calculating the phase angle as a function of 
spatial position of the first detector relative to any other detector and identifying the 
component at each detector having the calculated phase angle. The sets are determined 
for all of the voxels. A NMR image is obtained from the sets wherein each set of Fourier 
components comprises a three-dimensional map of the intensity variation of the NMR 

25 signal produced by each voxel of the magnetized body. This practice of the invention 
preferably employs a Fourier transform algorithm, described in Fourier Transform 
Reconstruction Algorithm Section, to determine the spatial location of each voxel from 
the corresponding set of Fourier components comprises a three-dimensional map of the 
intensity variation of the NMR signal produced by each voxel of the magnetized body 

30 over the sample space. This is repeated for each set to form the NMR image of the 
object. 

NMR images produced according to prior art methods and systems rely on 
applying an additional magnetic field in the direction of the primary field having a 
gradient along an axis in the transverse plane to cause a phase variation of the NMR 
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signal along the axis in the transverse plane. The direction axis of the gradient is varied 
a plurality of times to gives rise to an equivalent number of lines in the reconstructed 
image. In the present invention, the unique phase variation of the NMR signal is 
provided by the combination of L) the angle 6 suspended between the direction of the 
5 detector and the radial vector , the vector from the dipole to the detector, and 2.) the angle 
</> due to a separation distance r between a voxel and a detector given by the 
wavenumber of the RF field k times r . 

My prior inventions disclosed in U.S. Patent No. 5,073,858 and U.S. Patent 
Application No. 09/191,454, the complete disclosures of both which are incorporated 

10 herein by reference, are in part based on the realization that matter having a permeability 
different from that of free space distorts a magnetic flux applied thereto. This property is 
called magnetic susceptibility. An object, herein called a phantom, can be considered as 
a collection of small volume elements, herein referred to as voxels. When a magnetic 
field is applied to the phantom, each voxel generates a secondary magnetic field at the 

15 position of the voxel as well as external to the phantom. The strength of the secondary 
magnetic field varies according to the strength of the applied field, the magnetic 
susceptibility of the material within the voxel, and the distance of the external location 
relative to the voxel. For example, my U.S. Patent No. 5,073,858 and U.S. application 
No. 09/191,454 teache that the net magnetic flux at a point extrinsic to a phantom to 

20 which a magnetic field is applied, is a sum of the applied field and the external 

contributions from each of the voxels. The '858 patent further teaches sampling the 
external flux point by point and employing a reconstruction algorithm, to obtain the 
magnetic susceptibility of each voxel from the sampled external flux. 

Unlike the '858 patent that relies on a static response from a magnetized body to 

25 determine the magnetic susceptibility of the body, my invention disclosed in U.S. Patent 
Application No. 09/191,454 elicits a radiative response from a magnetized body by 
subjecting the body to a resonant radiation field. One embodiment of the '454 invention 
generates a three-dimensional magnetic susceptibility image of an object including a 
patient placed in a magnetic field from a three-dimensional map of a radio frequency 

30 (RF) magnetic field external to the patient, induced by subjecting selected nuclei of the 
body to a resonant RF field. Application of an RF pulse to the body causes the body to 
emit the RF magnetic flux external to the body. A Fourier transform of this external flux 
produces its frequency components ("Larmor frequencies"). Each Larmor frequency is 
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used to determine the magnetic susceptibility of the voxels of the body producing that 
Larmor frequency. Further, the intensity variation of the external RF field over a three- 
dimensional volume of space is used to determine the coordinate location of each voxel. 
One practice of the inventions disclosed in my U.S. Patent No. 5,073,858 and 
5 U.S. Patent Application No. 09/191,454 obtains a three-dimensional magnetic 

susceptibility map of a magnetized body from a three-dimensional map of a secondary 
magnetic flux produced by the magnetized body, and detected over a three-dimensional 
volume of space external to the body, herein referred to as the "sample space." The 
extrinsic magnetic flux is sampled at least at the Nyquist rate, i.e., at twice the spatial 

10 frequency of the highest frequency of the Fourier transform of the magnetic 

susceptibility map of the phantom, to allow adequate sampling of spatial variations of the 
external magnetic flux. This practice of the inventions preferably employs a Fourier 
transform algorithm, described in Fourier Transform Reconstruction Algorithm Section, 
to form the magnetic susceptibility map of the object. 

15 The present invention relates to systems for providing images of distributions of a 

quantity, in a chosen region of the body, by gyromagnetic resonance, particularly nuclear 
magnetic resonance (NMR) techniques. Such techniques may be used for examining 
bodies of different kinds. A particularly beneficial application is the examination of 
patients for medical purposes. Unlike my U.S. Patent No. 5,073,858 and U.S. Patent 

20 Application No. 09/191,454, the present invention employs nuclear magnetic resonance 
(NMR) to induce a magnetized phantom of essentially constant magnetic susceptibility 
to emit an external radiation having a magnetic field component. In particular, 
application of an RF pulse, resonant with selected nuclei of a magnetized body, can 
polarize the nuclei through rotation of their magnetic moments. The polarized nuclei 

25 within a voxel precess about the local magnetic field in the voxel at a Larmor frequency 
determined by the applied magnetic field at position of the voxel. The superposition of 
external RF fields produced by all the voxels of the body creates the total external RF 
field at each detector that is time dependent. The external RF field recorded at the 
detectors as a function of time contains components each having a unique phase angle 

30 relative to other components. Each component corresponds to an emitting voxel of the 
phantom. The time dependent signal at each detector maybe transformed into a series of 
components having intensity and phase data. Each set of components of the NMR signal 
over the sample space due to a given voxel is determined from the phase data and the 
detector positions. The spatial variation of the NMR signal over the sample space is 
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used to determine the location of the voxel in the phantom. This is repeated for all sets 
of components, each corresponding to a voxel to reconstruct the NMR image. 

The radiation source for magnetizing a body to be imaged can be a direct current 
("DC") magnet, including a superconducting magnet. The radiation sources and 
5 amplifiers for applying an RF pulse to the magnetized body are well known in the art, 
and include, but are not limited to, klystrons, backward wave oscillators, Gunn diodes, 
Traveling Wave Tube amplifiers. A preferred embodiment of the invention employs a 
three dimensional array of antennas as detectors for sensing the external RF field. 

One practice of the invention detects the external RF field in the near field region 

10 where the distance of a detector sensing radiation from a voxel at a distance r from the 
detector is much smaller than the wavelength X of the radiation emitted by the voxel, 
i.e., r « X (or kr « 1). The near fields are quasi-stationary, that is they oscillate 
harmonically as e~ mt , but are otherwise static in character. Thus, the transverse RF 
magnetic field of each voxel is that of a dipole. In one embodiment, an array of 

15 miniature RF antennas sample the external RF field over a three-dimensional volume of 
space that can be either above or below the object to be imaged. The distance r of the 
detector from the voxel gives rise to the phase term e~ lkr of the component of the 
detected RF signal where k is the wavenumber of the NMR signal. The harmonic 
oscillation of each RF dipole is equivalent to the dipole rotating in the transverse plane. 

20 The detector is responsive to a component in this plane. At a point in time, each RF 

dipole is directed at an angle 0 relative to the direction of detection of the detector. The 
phase angle 6 of the RF dipole relative to the direction of detection axis of the detector 
gives rise to a phase angle term e~ l9 . In a preferred embodiment, the sum of the phase 
angles, kr and 9 , are unique for each voxel at each detector. The position of each 

25 detector relative to a different detector may be used to calculated the phase angle of the 
second relative to the first. This may be repeated over all of the detectors to give the set 
of intensities of the NMR signal over the sample space due to a voxel. The location of 
each voxel is determined through the spatial variations of the intensity of the NMR field 
of the set of components associated by phase. Thus, the phase of the components of the 

30 external RF radiation, and the intensity variations of the external RF radiation provide 
the necessary information for providing a NMR image of the magnetized phantom, such 
as a human body. Such a NMR image can be employed to obtain anatomical images of a 
human body based on selected physiological parameters. 
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Li an embodiment of the present invention, the NMR image of an object 
including a patient placed in a magnetic field is generated from a three-dimensional map 
of the transverse resonant radio frequency (RF) magnetic flux external to the patient. 
The external RF field recorded at the detectors as a function of time contains components 
5 each having a unique phase angle relative to other components. Each component 
corresponds to an emitting voxel of the phantom. The time dependent signal at each 
detector may be transformed into a series of components having intensity and phase data. 
Each set of components of the NMR signal over the sample space due to a given voxel is 
determined from the phase data and the detector positions. The intensity variation of the 
10 transverse RF field over the sample space is used to determine the coordinate location of 
each voxel. The RF field is the near field which is a dipole that serves as a basis element 
to form a unique reconstruction. The geometric system function corresponding to a 
dipole which determines the spatial intensity variations of the RF field is a band-pass for 
k p =k z . Preferably, each volume element is reconstructed independently in parallel with 

15 all other volume elements such that the scan time is no greater than the nuclear free 
induction decay (FID) time. 

Secondary Magnetic Field 

20 The magnetic moment m z of each voxel is a magnetic dipole. And the phantom 

can be considered to be a three-dimensional array of magnetic dipoles. At any point 
extrinsic to the phantom, the z-component of the secondary flux, B , from any single 
voxel is 

2z 2 -x 2 ~y 2 
(x +y +z ) 

25 where x, y, and z are the distances from the center of the voxel to the sampling point. It 
is shown in APPENDICES I-IV that no geometric distribution of magnetic dipoles can 
give rise to Eq. (1). Therefore, the flux of each magnetic dipole (voxel contribution) 
forms a basis set for the flux of the array of dipoles which comprise the NMR image of 
the phantom. 

30 Eq. (1) is a system function which gives the magnetic flux output in response to a 

magnetic dipole input at the origin. The phantom is an array of spatially advanced and 
delayed dipoles weighted according to the magnetic moment of each voxel; this is the 
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input function. The secondary flux is the superposition of spatially advanced and 
delayed flux, according to Eq. (1); this is the output function. Thus, the response of 
space to a magnetized phantom is given by the convolution of Eq. (1) with the series of 
weighted, spatially advanced and delayed dipoles representing the NMR image of the 
5 phantom. 

In Fourier space, the output function is the product of the Fourier transform (FT) 
of the system function and the FT of the input function. Thus, the system function filters 
the input function. The output function is the flux over all space. However, virtually all 
of the spectrum (information needed to reconstruct the NMR image) of the phantom 
10 exists in the space outside of the phantom because the system function is essentially a 
band-pass filter. This can be appreciated by considering the FT, H\k p ?k z ~\, of Eq. (1): 

where k p is the spatial frequency in the xy-plane or k p -plane and k z is the spatial 
frequency along the z-axis. H[k p ,k 2 ] is a constant for k p and k z essentially equal as 
15 demonstrated graphically in FIGURE lc. 



Band-Pass Filter 

When a static magnetic field H 0 with lines in the direction of the z-axis is applied 
20 to an object comprising a material containing nuclei such as protons that possess 
magnetic moments, the field magnetizes the material. As a result a secondary field 
superposes the applied field as shown in FIGURE 9. In the applied magnetic field, the 
magnetic moments of each nuclei precesses about the applied magnetic field. However, 
the magnetization of any one nucleus is not observed from the macroscopic sample. 
25 Rather the vector sum of the dipole moments from all magnetic nuclei in the sample is 
observed. This bulk magnetization is denoted by the vector M . In thermal equilibrium 
with the primary field H 0 , the bulk magnetization M is parallel to H 0 . The 
magnetization vector then comprises magnetic dipole m . The secondary field outside of 
the object (phantom) and detected at a detector 301 is that of a series of magnetic dipoles 
30 centered on volume elements 302 of the magnetized material. In Cartesian coordinates, 
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the secondary magnetic flux, B 9 at the point (x,y,z) due to a magnetic dipole having a 
magnetic dipole moment m z at the position (x 09 y 09 z 0 )' is 

_ m s (2(z - zof - (x - xo) 2 - (y -y 0 ) 2 \ 
[(x-xo) ±(y-y 0 ) +(2-zo)r 

(2z 2 -x 2 -y 2 } 

5 B ^ r V 2, 2 , 2l5 ^ ®^5(x-x 0? ;;-y 03 z-z 0 X ( 4 ) 

ix -t y + z J 

where i z is the unit vector along the z-axis. The field is the convolution of the system 
function, h(x,y 9 z) or h(p 9 <j) 9 z) (the left-handed part of Eq. (4)), with the delta function 
(the right-hand part of Eq. (4)), at the position (x 0 9 y 09 z 0 ) . A very important theorem of 
Fourier analysis states that the Fourier transform of a convolution is the product of the 
10 individual Fourier transforms [2]. The Fourier transform of the system function, 
h(x 9 y,z) orh(p 9 ^,z), is given in APPENDIX V. 

The z-component of a magnetic dipole oriented in the z-direction has the system 
function, h(x 9 y 9 z) 9 which has the Fourier transform, H\k x9 k 9 k z ] 9 which is shown in 
FIGURE lc. 

'Tit (6) 

V 

The output function, the secondary magnetic field, is the convolution of the 
system function, h(x 9 y 9 z)~~ the geometric transfer function for the z-component of a z- 
oriented magnetic dipole with the input function— a periodic array of delta functions each 
20 at the position of a magnetic dipole corresponding to a magnetized volume element. 
(2z 2 -x 2 ~y 2 } 

,2 2 2~JH® ^tn z 8(x-nx Q ,y-ny 0 ,z-nz<>) (7) 

[x +y +z J 

The Fourier transform of a periodic array of delta functions (the right-hand side of Eq. 
(7)) is also a periodic array of delta functions in k-space: 

— t^Jk x ~-,k y --,K--) (8) 



Xoy Q Z Q n=-cc 
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By the Fourier Theorem, the Fourier transform of the spatial output function, Eq. (7), is 
the product of the Fourier transform of the system function given by Eq. (6), and the 
Fourier transform of the input function given by Eq. (8). 
4k 1 ^ J , n , n , n] 



5 In the special case that 

K=K (io) 

the Fourier transform of the system function (the left-hand side of Eq. (9)) is given by 

H=4ti (11) 
Thus, the Fourier transform of the system function band-passes the Fourier transform of 

10 the input function. Both the input function (the right-hand part of Eq. (7)) and its Fourier 
transform (the right-hand part of Eq. (9)) are a periodic array of delta functions. No 
frequencies of the Fourier transform of the input function are attenuated; thus, no 
information is lost in the case where Eq. (10) holds. 

In an embodiment of the present invention, the magnetization vector is rotated 

15 into the transverse plane by an additional RF field Hj . The magnetization vector then 
comprises a rotating magnetic dipole m in the transverse plane. The NMR image may 
be reconstructed by sampling the external field from a series of RF dipoles rather than 
that from a series of static dipoles. In this case, the Fourier transform of the system 
function also band-passes the Fourier transform of the input function. Thus, the 

20 resolution of the reconstructed NMR image is limited by the spatial sampling rate of the 
secondary RF magnetic field according to the Nyquist Sampling Theorem. 



Reconstruction 

25 The NMR image may be reconstituted using a Fourier transform algorithm. The 

algorithm is based on a closed-form solution of the inverse problem — solving the spatial 
distribution of an array of magnetic dipoles from the measured extrinsic secondary (RF) 
field that is transverse to the magnetic flux that magnetizes the voxels. The transverse 
RF magnetic field of each voxel is that of a dipole, the maximum amplitude is given by 

30 Eq. (1) wherein the Larmor frequency of each voxel is essentially the same, and m z , the 
magnetic moment along the z-axis, of Eq. (1) corresponds to the bulk magnetization M 
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of each voxel. In terms of the coordinates of Eq. (1), an array of miniature RF antennas 
point samples the maximum dipole component of the RF signal over the sample space 
such as the half space above (below) the object to be imaged wherein each RF signal is 
has a unique phase shift due to the relative spatial relationship of each voxel and each 
5 detector. The phase shift of each component and the relative spatial relationship of the 
detectors is used to assign a component from each detector to a set. Each set of 
components associated by phase comprises the spatial variations of the intensity of the 
transverse RF field of any given voxel. The intensity variation over the sample space is 
used to determine the coordinate location of each voxel. In the limit, each volume 
10 element is reconstructed independently in parallel with all other volume elements such 
that the scan time is no greater than the nuclear free induction decay (FID) time. 

The NMR scan performed on the object to be imaged including a human 
comprises the following steps: 

15 • The magnetic moments of nuclei including protons of the object to be imaged that are 
aligned by the primary field are further aligned by a radio frequency (RF) pulse or 
series of pulses. 

• The strength and duration of the rotating H l (RF) field that is resonant with the protons 
20 of the magnetized volume and is oriented perpendicularly to the direction of the 

magnetizing field is applied such that the final precession angle of the magnetization 
is 90 ° (0 Hi = 90° ) such that the RF dipole is transverse to the primary magnetizing 

field and perpendicular to the RF magnetic field detector. 

25 • NMR pulse sequences which provide the signals for a T x or T 2 image may be applied. 

For example a 90 ° pulse may be followed by a series of 180 ° pulses. One sequence 
is the Carr-Purcell-Meiboom-Gill (CPMG) sequence [3]. 

• The free induction decay signals are recorded. 

30 

• The time dependent signals are Fourier transformed to give the intensity and phase of 

each component. The NMR signal of each voxel at any given detector gives rise to a 
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Fourier component with a unique phase angle relative to the Fourier component of 
any other voxel of the phantom at that detector. 



• The matrix of Fourier components that correspond to the NMR signal of a given voxel 
5 over the detectors is determined. This may be achieved by using a first component 
having a phase angle and calculating the phase angle as a function of spatial position 
of the first detector relative to any other detector and identifying the component at 
each detector having the calculated phase angle. The matrices are determined for all 
of the voxels. The measurements of the spatial variations of the transverse RF field 
10 of a given matrix is used to determine the coordinate location of each voxel. Thus, 

each matrix of components associated by phase comprises the intensity variation over 
the sample space of the RF field of the bulk magnetization M of each voxel. 



• The Fourier transform algorithm is performed on each set of components over the 
detector array to map each bulk magnetization M corresponding to a voxel to a 
spatial location over the image space. The bulk magnetization map (NMR image; 
.also the input function) is given as follows. With respect to the coordinate system of 
Eq. (1), (x, y, and z are the Cartesian coordinates, rn z , the magnetic moment along 
the z-axis, of Eq. (1) corresponds to the bulk magnetization M of each voxel, and B 
is the magnetic flux due to the magnetic moment shown in FIGURE 9; the 
relationship to the NMR coordinate system is given in the Reconstruction Algorithm 
Section) the origin of the coordinate system, (0,0,0), is the center of the upper edge of 
the phantom. The phantom occupies the space below the plane x, y, z = 0 (z < 0 in 
the phantom space), and the sampling points lie above the plane ( z > 0 in the 
sampling space). The magnetic flux in the sampling space is given by multiplying 
the convolution of the input function with the system function by the unitary function 
(one for z > 0 and zero for z < 0). The input function can be solved in closed-form 
via the following operations: 

1 . Record the RF NMR signal at discrete points in the sampling space. Each point is 
designated (x, y, z, RF) and each RF value is an element in matrix A. The time 
dependent signals are Fourier transformed to give the intensity and phase of each 
component. The NMR signal of each voxel at any given detector gives rise to a Fourier 
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component with a unique phase angle relative to the Fourier component of any other 
voxel of the phantom at that detector. The matrix of Fourier components that correspond 
to the NMR signal of a given voxel over the detectors is determined. This may be 
achieved by using a first component having a phase angle and calculating the phase angle 
5 as a function of spatial position of the first detector relative to any other detector and 
identifying the component at each detector having the calculated phase angle. The 
matrices A n are determined for all of the voxels. The measurements of the spatial 
variations of the transverse RF field of a given matrix is used to determine the coordinate 
location of each voxel. Thus, each matrix of components associated by phase comprises 
10 the intensity variation over the sample space of the RF field of the bulk magnetization M 
of each voxel. 

2. Discrete Fourier transform each matrix A n to obtain each matrix B n . 

15 3. Multiply each element of each matrix B n by the corresponding inverse (reciprocal) 
value of the Fourier transform of the system function, Eq. (1), evaluated at the same 
frequency as the element of the matrix A n . This is matrix C n . 

4. Generate matrix D n by taking the discrete inverse Fourier transform of matrix C n . 

20 

5. Multiply each element of each matrix D n by the distance squared along the z-axis to 
which the element corresponds to generate the position of the bulk magnetization M of 
voxel n. (This corrects the limitation of the sample space to z > 0). 

25 6. In one embodiment, the point spread of the reconstructed voxel is corrected by 

assigning one voxel above a certain threshold with the bulk magnetization M . The other 
voxels are assigned a zero value. 

7. This procedure is repeated for all matrices A n . In the limit with sufficient phase 
30 resolution, each volume element is reconstructed independently in parallel with all other 
volume elements such that the scan time is no greater than the nuclear free induction 
decay (FID) time. 
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BRIEF DESCRIPTION OF THE DRAWINGS 



FIGURE la is a schematic diagram of an apparatus according to the teachings of 
the invention in which a magnet magnetizes a body to be imaged, an RF generator 
5 excites the body, and an array of detectors detect radiation emitted by the body in 
response to the RF excitation radiation; 

FIGURE lb is a schematic diagram of an apparatus according to the teachings of 
the invention that employs a three-dimensional array of detectors to detect radiation 
emitted by nuclei of a body to be imaged in response to excitation of the nuclei by a 
10 radiation source; 

FIGURE lc is a plot of the Fourier transform H[k x9 k y9 k z ] of the system function 

h(x,y 9 z) (Eq. (1)) corresponding to the z-component of a magnetic dipole oriented in the 
z-direction in accordance with the invention; 

FIGURE 2 is the plot of the field of a ring of dipoles of radius R and magnetic 
15 moment m = 10 4 Gem 3 given by Eq. (1.14) as a function of radius R where the position 
of the center of the ring relative to the detector is the point (0,0,10) in accordance with 
the invention; 

FIGURE 3 is the plot of the field of a ring of dipoles of radius R = 0.2 cm and 
magnetic moment m = 10 4 Gem given by Eq. (1.14) as a function of the distance 
20 between the detector at the origin and the center of the ring at the points (0,0, 
z = 4 cm to z = 1 5 em ) in accordance with the invention; 

FIGURE 4 is the plot of the field of a shell of dipoles of radius R and magnetic 
moment m = 10 4 Gem given by Eq. (11.17) as a function of radius R where the position 
of the center of the shell relative to the detector is the point (0,0,10) in accordance with 
25 the invention; 

FIGURE 5 is the plot of the field of a shell of dipoles of radius R = 0.2 cm and 
magnetic moment m = 10 4 Gem 3 given by Eq. (11.17) as a function of the distance 
between the detector at the origin and the center of the shell at the points (0,0, 
z = 4 cm to z = 15 cm ) in accordance with the invention; 
30 FIGURE 6 is the plot of the field of a sphere of dipoles of radius R and magnetic 

moment m = 1 0 4 Gem 3 given by Eq. (IV. 16) as a function of radius R where the 
position of the center of the sphere relative to the detector is the point (0,0,10) in 
accordance with the invention; 
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FIGURE 7 is the plot of the field of a sphere of dipoles of radius R = 0.2 cm and 
magnetic moment m = 10 4 Gem given by Eq. (IV. 16) as a function of the distance 
between the detector at the origin and the center of the sphere at the points (0,0, 
z = 4 cm to z = 15 cm ) in accordance with the invention; 
5 FIGURE 8 shows a typical the nuclear magnetic resonance (4D-MRI) apparatus 

in accordance with the invention; 

FIGURE 9 shows the coordinate system (x, y, z) of Eq. (1) with a primary field 
H 0 and the corresponding magnetic dipole both oriented parallel to the z-axis wherein 
the z-component of the flux due to the z-oriented dipole is measured at a detector 
10 according to Eq. (1) and shows the distances from the voxel to the detector in accordance 
with the invention; 

FIGURE 10 shows the stationary coordinate system of the nuclear magnetic 
resonance (4D-MRI) apparatus of FIGURE 8 corresponding to the coordinate system of 
FIGURE 9 in accordance with the invention; 

1 5 FIGURE 1 1 shows the rotating NMR coordinate system (x R 5 y R , z R ) and the 

stationary coordinate system (x, y, z) of the NMR detector corresponding to the 
coordinate system of FIGURE 9 and FIGURE 10 of a primary field H 0 oriented parallel 
to the z R -axis and the z-axis and the corresponding transverse RF magnetic dipole 
oriented in the x R y R -plane and periodically parallel to the y-axis wherein the spatial 

20 variation of the RF y-component of the flux due to the RF dipole is measured at a 

detector according to Eq. (1) and shows the distances from the voxel to the detector in 
accordance with the invention; 

FIGURE 12 is a schematic of the three dimensional detector (antennae) array 
with respect to the stationary NMR coordinate system of FIGURE 10 which corresponds 

25 to the coordinate systems shown in FIGURES 9 and 1 1 in accordance with the invention; 
FIGURE 13 shows the general process of reconstruction by reiteration, and 
FIGURE 14 shows the stationary coordinate system (x, y 5 z) of the NMR detector 
corresponding to the coordinate system of FIGURE 9 and FIGURE 10 of a primary field 
H 0 oriented parallel to the z R -axis and the z-axis and the corresponding transverse RF 

30 magnetic dipole oriented in the x R y R -plane and periodically parallel to the y-axis 

wherein the spatial variation of the RF y-component of the flux due to the RF dipole is 
measured at a detector according to Eq. (1) and shows the distances and angles between a 



WO 02/16956 PCT/US01/25954 

16 

voxel linear to a first detector., a second nonaligned voxel, and a second detector in 
accordance with the invention; 

Further details regarding specific derivations, and calculations are provided in the 
attached appendices, wherein: 
5 APPENDIX I is the field produced by a ring of dipoles according to the present 

invention; 

APPENDIX II is the derivation of the field produced by a shell of dipoles 
according to the present invention; 

APPENDIX III is the mathematical proof that the field produced by a shell of 
10 magnetic dipoles is different from that of a single dipole according to the present 
invention; 

APPENDIX IV is the derivation of the field produced by a sphere of dipoles 
according to the present invention; 

APPENDIX V is the derivation of the Fourier transform of the system function of 
15 the z-component of the magnetic field from a dipole oriented in the direction of the z- 
axis used in the reconstruction process according to the present invention; 

APPENDIX VI is the derivation of S = HF® U(k z ) convolution from Eq. (55) 
used in a reconstruction process according to the present invention, and 

APPENDIX VII is the derivation of the Inverse Transform of Eq. (69) to Give 
20 Inverse Transform 1, Eq. (69), used in a reconstruction process according to the present 
invention. 

DETAILED DESCRIPTION OF THE INVENTION 

25 An exemplary embodiment of a nuclear magnetic resonance (4D-MRI) apparatus 

110 according to the teachings of the present invention is shown in FIGURE la. The 
apparatus 110 includes a magnet 112, such as a superconducting magnet, that provides a 
primary or magnetizing field, to magnetize a body 1 14 to be imaged. A magnetometer 
116 determines the primary or magnetizing field in the volume to be occupied by the 

30 body, i.e., in the image space, in the absence of the body. One practice of the invention 
utilizes a magnetometer that employs NMR of protons in water for determining the 
primary radiation field at multiple points in the image space. In an embodiment, the 
primary field is uniform as recorded by magnetometer 116. 
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A radiation source 1 18, such as a radio frequency generator, applies an RF pulse 
to the body in combination with RF antennae 1 18a when the body is placed in the 
magnetic field, to excite and thereby polarize selected nuclei of the body. The excited 
nuclei emit an RF radiation that a plurality of detectors 120, disposed in a plane above or 
5 below the object, can detect. The excitation pulse can be selected to rotate the 

magnetization of the nuclei, preferably by 90 degrees, with the respect to the primary 
field. In such a case, the RF radiation that the excited nuclei emit is primarily along a 
direction perpendicular to the plane of the detectors. In an embodiment using techniques 
known by those skilled in the art, RF pulse sequences are applied to generate the data for 

10 7\ and/or T 2 NMR images. The detectors can be selected to respond only to components 
of a magnetic field perpendicular to the plane in which they reside. Thus, such detectors 
can detect the emitted RF field without interference from other components of the 
magnetic field to permit a unique reconstruction of the NMR image. In an embodiment, 
a drive mechanism 122 moves the detectors 120 in a direction perpendicular to the plane 

15 of the detectors, to sample the external RF field over a three-dimensional volume. The 
separation of the detectors and the step size of the movement of the detectors along a 
direction perpendicular to the plane are selected such that the detectors sample the 
external RF field over a three-dimensional volume, i.e., the sample volume, at least at the 
Nyquist rate. Preferably, a three dimension array of detectors is used to sample the RF 

20 field at the Nyquist rate. Such a sampling advantageously allows obtaining the NMR 
image of the body. An embodiment, employs an impedance-matched array of RF 
antennas that are time-multiplexed to reduce cross talk among them. The RF field may 
be sampled synchronously or the field may be sampled at known times so that the phase 
at any given detector may be related to that at any other detector. 

25 The process of magnetization of selected nuclei of the body can be better 

understood by referring to FIGURE 1 1 that shows a magnetizing field H 0 applied to a 
voxel 14a of the body 1 14. The RF excitation field is selected to be in a direction 
perpendicular to the field H 0 , and its magnitude is designated as Hj . NMR active nuclei 
of the voxel, such as protons, possess both angular momentum and a magnetic moment. 

30 Thus, the vector sum of the magnetic moments of all such NMR active nuclei present in 
the voxel 14a give rise to the bulk magnetization of the voxel 14a. The bulk 
magnetization vector M (not shown) designates the collective contribution of a selected 
type of NMR active nuclei to the magnetization of the voxel 14a which corresponds to a 
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RF parallel magnetic moment ra . The RF excitation field is selected to be in resonance 
with the selected type of nuclei, e.g., protons, of the voxel 14a, to rotate the bulk 
magnetization vector M . In a rotating frame representation, designated in FIGURE 1 1 
as ( x r ■> y R . z r) and well known to those skilled in NMR, the magnetization vector M 
5 rotates about H 1 so long as the RF radiation is present. The rate of rotation of the 
magnetization vector about the applied field H t depends on the gyromagnetic ratio of 
the affected nuclei and the magnitude of the excitation field, U t . The duration of the RF 
radiation can thus be selected to cause a rotation of the magnetization vector, initially 
aligned along H 0 , onto the x R y R -plane, for example by a 90° rotation. After the RF 
10 excitation field is turned off, the rotated magnetization M precesses about H 0 at the 
position of the voxel 14a. 

The precession of the magnetization M of each voxel about H 0 produces a 
radiating dipolar field corresponding to a magnetic moment m that has its maximum 
intensity along the y R direction. The precession frequency of the magnetic moment 
1 5 called the Larmor frequency is dependent on the magnetic flux of the primary field 

which is uniform in a preferred embodiment. A detector 20a can be, for example, an RF 
antenna pointing along the y direction to respond selectively to the radiating RF field 
from all of the voxels. This external RF field is recorded as a function of time at each 
over the entire sample space. The NMR signal at each detector as a function of time is 
20 processed by a method such as a Fourier transform operation to give a plurality of 
Fourier components each having an intensity and a phase angle. The NMR signal of 
each voxel at any given detector gives rise to a Fourier component with a unique phase 
angle relative to the Fourier component of any other voxel of the phantom at that 
detector. Each set of components of the NMR signal over the sample space due to a 
25 given voxel is determined from the phase data and the detector positions. The spatial 
variation of the NMR signal over the sample space is used to determine the location of 
the voxel in the phantom. This is repeated for all sets of components, each 
corresponding to a voxel to reconstruct the NMR image. Preferably, the position of each 
voxel is reconstructed independently in parallel with all other voxels such that the scan 
30 time is no greater the time for the excited nuclei to return to their unexcited state called 
the nuclear free induction decay (FID) time. This may be achieved by using sufficiently 
high spatial and temporal sampling rate such that each set of components associated by 
phase corresponds to a single voxel. 
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With further reference to FIGURE la, an analog to digital converter ("A/D 
converter") 124 converts the analog signal outputs of the detectors 120 into digital 
signals. A preferred embodiment of the invention utilizes at least a 12 bit A/D converter 
to digitize the output signals. A reconstruction processor 126 receives the digital signals, 
5 and determines the NMR image of the body. The reconstruction processor 126 includes 
a Fourier transform processor 126a that obtains the phase components of the external RF 
field. In addition, a spatial locator 126b, which is a part of the reconstruction processor 
126, employs the variation of the maximum intensity of the external RF field over the 
sample space (the three dimensional space sampled by the detectors) to locate the 

10 positions of the voxels in the image space producing a particular set of components of 
the external RF field associated by phase, in a manner described in detail below (See the 
Reconstruction Algorithm Section). An display processor 128 displays a two- 
dimensional or a three-dimensional image corresponding to the NMR image of the body. 
The spatial locator 126b employs an algorithm, described in detail below in the 

15 Reconstruction Algorithm Section, to determine the positions of the voxels of the 

magnetized body that produce each set of components of the external RF field associated 
by phase. It is demonstrated in the Uniqueness of Reconstruction Section that a RF 
magnetic field produced by a geometric distribution of dipoles is unique. Therefore, a 
unique spatial distribution of magnetic dipoles, such as those corresponding to the bulk 

20 magnetization of each voxel due to precessing nuclei, gives rise to a unique magnetic 
field. Thus, the measured external RF field can provide a unique solution for the spatial 
distribution of magnetic dipoles, i.e., magnetized voxels comprising excited nuclei, in the 
body. 

Referring again to FIGURE la, the spatial locator preferably employs a Fourier 
25 Transform Reconstruction Algorithm, described in detail in the Reconstruction 

Algorithm Section, to map a bulk magnetic moment M corresponding to a particular set 
of components associated by phase of the external RF field onto a spatial location or 
locations in the body. The latter case applies if more than one location in the body gives 
rise to a particular phase component of the external RF radiation. 
30 As shown in FIGURE la, a preferred embodiment of the invention employs an 

open design magnet, such as a Helmholtz coil design, to allow positioning the array of 
detectors 120 close to the body 1 14 to be imaged. The orientation of the magnetic field 
with respect to the body can be selected to optimize the signal to noise ratio of the 
signals detected by the array of detectors 120. For example, in case of imaging a patient 
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body, the primary magnetic field can be selected to be coaxial with the body, or it can 
alternatively be perpendicular to the body axis. 

The 4D-MRI apparatus of FIGURE la provides a number of advantages. In 
particular, because all of the data is acquired at once rather than over hundreds of 
5 repeated sequences of pulses as is the case with NMR systems of the prior art, the 
apparatus of tjie invention is particularly suited for imaging cardiopulmonary and 
vascular systems. Further, since the apparatus may increase the number of detectors to 
increase the resolution, the apparatus of the invention may achieve a higher resolution 
such as 10" 3 cm 3 ; thereby, permitting physicians to view human anatomy and pathology 
10 in a manner not available with conventional imaging techniques. Further, the present 
technique may provide a three dimensional image that can be displayed from any 
perspective. 

An alternative embodiment of the nuclear magnetic resonance (4D-MR1) 
apparatus of the invention, shown in FIGURE lb, employs a three-dimensional array of 

15 detectors 230, spaced apart to detect spatial variations of the emitted RF radiation at least 
at the Nyquist frequency. A magnet 212 provides a magnetizing in a volume to be 
occupied by a body, i.e., image space. A magnetometer 216 measures the magnetizing 
field at a plurality of positions in the image space in the absence of the body to determine 
the uniformity of the primary field. As in embodiment of FIGURE la, the body 214 to 

20 be imaged is placed in a magnetizing field provided by the magnet 212. An RF 

generator 21 8 in combination with an RF antenna 218a apply an RF pulse or a sequence 
of RF pulses to the body to polarize selected nuclei of the body. The three-dimensional 
array of detectors 230 provide output signals in response to the RF radiation emitted by 
the body. The RF signal may be recorded synchronously to permit the relative phases of 

25 the Fourier components comprising the RF signal. In an embodiment, the detection may 
be synchronized relative to the excitation. A digitizer 224 digitizes the output signals 
and sends the digital signals to a construction processor 226 that determines variations of 
the NMR image of the body in a manner similar to that described in connection with the 
embodiment of FIGURE la. A display processor 228 receives the information regarding 

30 the spatial variations of the intensity of the NMR signal from the construction processor 
226, and provides a two-dimensional or a three-dimensional NMR image of the body. 
Employing a three-dimensional array of detectors advantageously decreases the 
acquisition time because the emitted RF signal over the entire three dimensional sample 
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space is detected at once. The shortening of the acquisition time in turn reduces motion 
artifacts in the NMR image. 

Uniqueness of Reconstruction 

5 

The nature of the RF field can be determined from Maxwell's equations applied 
to a sinusoidal current. With a sinusoidal current J(x f ) confined to small region 
compared with a wavelength, the solution of the vector potential A(x) is [4] 

A(x) = ^J(xO^^ (12) 

10 where k = — is the wavenumber, and a sinusoidal time dependence is understood. The 

c 

magnetic induction is given by 

B = VxA (13) 
while, outside the source, the electric field is 

E = -VxB (14) 

k 

For a source of dimension d , the fields in the far zone defined by d « X « r are 
transverse to the radius vector and fall off as r~\ typical of radiation fields. For the near 
zone where r « X (or kr « 1 ), the exponential in Eq. (12) can be replaced by unity. 
Then the vector potential is given by 

/cr — > 0 c lm 2/ + 1 r J 

This shows that the near fields are quasi-stationary, oscillating harmonically as e~ mt 9 but 
otherwise static in character. 

Nuclear magnetic resonance (NMR), which is commonly called magnetic 
resonance imaging (MRI), is a means to measure the primary and secondary magnetic 
fields to provide the input to the NMR reconstruction algorithm. The proton 
gyromagnetic ratio y p I 2n is 

y p I 2tv =42.57602 MHz T~ l (16) 

The NMR frequency / is the product of the proton gyromagnetic ratio given by Eq. (16) 
and the magnetic flux B . 

/ = r P I 2nB= 42. 57602 MHz r l B (17) 
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A typical flux for a superconducting NMR imaging magnet is 0. 25 T . According to Eq. 
(17) this corresponds to a radio frequency (RF) of 10.6 MHz which corresponds to a 
wavelength of 28.3 m . In the present invention, each RF antennae of an array is located 
at a distance of about 10 cm from the voxels within the image space. Thus, the RF field 
5 is detected in the near zone where r « X (or kr « 1 ), and the near fields according to 
Eq. (15) are quasi-stationary, oscillating harmonically as e mt , but otherwise static in 
character. The transverse RF magnetic field of each voxel is that of a RF dipole, the 
maximum amplitude is given by Eq. (1) wherein the Larmor frequency of each voxel is 
determine by the uniform primary field H 0 , and m z , the magnetic moment along the z- 

10 axis, of Eq. (1) corresponds to the bulk magnetization M of each voxel. 

An object containing nuclei with a magnetic moment, herein called a phantom, 
can be considered as a collection of small volume elements or voxels. When a static 
magnetic field H 0 with lines in the direction of the z-axis is applied to an object 
comprising a material containing nuclei such as protons that possess magnetic moments, 

15 the field magnetizes the material. As a result a secondary field superposes the applied 
field as shown in FIGURE 9. In the applied magnetic field, the magnetic moments of 
each nuclei precesses about the applied magnetic field. However, the magnetization of 
any one nucleus is not observed from the macroscopic sample. Rather the vector sum of 
the dipole moments from all magnetic nuclei in the sample is observed. This bulk 

20 magnetization is denoted by the vector M . In thermal equilibrium with the primary field 
H 0 , the bulk magnetization M is parallel to H 0 . In an embodiment of the present 
invention, the magnetization vector is rotated into the transverse plane by an additional 
RF field Hj . The magnetization vector then comprises a rotating magnetic dipole m in 
the transverse plane. The NMR image may be reconstructed by sampling the external 

25 field from a series of RF dipoles. 

T he field strength of a magnetic dipole moment is a function of the external 
position in space relative to the dipole. For convenience of analysis, the field of a series 
of static dipoles m having the coordinates shown in FIGURE 9 is analyzed for 
uniqueness. (The uniqueness of the field of a set of static dipoles applies for the 

30 equivalent case of RF dipoles oriented in the transverse plane.) Considering FIGURE 9, 
the net magnetic field at a point extrinsic to the phantom is a sum of the applied field and 
the contributions of each of the voxels within the object, the secondary field. The field is 
point sampled over a three dimensional space. 
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The secondary magnetic field due to magnetized tissue has to be modeled as 
noninteracting dipoles aligned with the imposed field. It is demonstrated below that the 
field of any geometric distribution of dipoles is unique, and the superposition principle 
holds for magnetic fields; therefore, a unique spatial distribution of dipoles gives rise to a 
5 unique secondary magnetic field, and it is further demonstrated below that this secondary 
field can be used to solve for the NMR image in closed form. It follows that this map is 
a unique solution. To prove that any geometric distribution of dipoles has a unique field, 
it must be demonstrated that the field produced by a dipole can serve as a mathematical 
basis for any distribution of dipoles. This is equivalent to proving that no geometric 

10 distribution of dipoles can produce a field which is identical to the field of a dipole. 

By symmetry considerations, only three distributions of uniform dipoles need to 
be considered. A magnetic dipole has a field that is cylindrically symmetrical. A ring, a 
shell, a cylinder, and a sphere of dipoles are the only cases which have this symmetry. A 
cylinder is a linear combination of rings. Thus, the uniqueness of the dipole field is 

15 demonstrated by showing that it is different from that of a ring, a shell, and a sphere. 
The uniqueness of the dipole from the cases of a ring, a shell, and a sphere of dipoles is 
demonstrated in APPENDIX I, APPENDIX II, and Appendix IV, respectively. The plot 
of the three cases of the field of a ring, shell, and a sphere of dipoles each of radius R 
and magnetic moment m = 10 4 Gem 3 given by Eq. (1.14), Eq. (II. 17), and Eq. (IV. 16) of 

20 APPENDIX I, APPENDIX II, and APPENDIX IV as a function of radius R where the 
position of the center of each distribution relative to the detector is the point (0,0,10) is 
given in FIGURES 2, 4, and 6, respectively. Since the fields vary as a function of radius 
R , the dipole field is not equivalent to these distributions of dipoles. It is further 
mathematically proven in APPENDIX III that the field produced by a shell of magnetic 

25 dipoles is different from that of a single dipole. All other fields are a linear combination 
of dipoles. Thus, the dipole is a basis element for the reconstruction of a NMR image. 
Since each dipole to be mapped gives rise to a unique field and since the total field at a 
detector is the superposition of the individual unique dipole fields, linear independence is 
assured; therefore, the NMR image is unique. In other words, there is only one solution 

30 of the NMR image for a given set of detector values which spatially measure the 

superposition of the unique fields of the dipoles. This map can be reconstructed using 
the algorithms described in the Reconstruction Algorithm Section. 

Eq. (1) is a system function which gives the magnetic flux output in response to a 
magnetic dipole input at the origin. The phantom is an array of spatially advanced and 
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delayed dipoles weighted according to the bulk magnetization of each voxel; this is the 
input function. The secondary flux is the superposition of spatially advanced and 
delayed flux, according to Eq. (1); this is the output function. Thus, the response of 
space to a magnetized phantom is given by the convolution of Eq. (1) with the series of 
5 weighted, spatially advanced and delayed dipoles representing the bulk magnetization 
map or NMR image of the phantom. The discrete signals are recorded by a detector 
array over the sample space comprising the xy-plane and the positive z-axis of FIGURE 
9. In an embodiment of the present invention, the magnetization vector is rotated into 
the transverse plane by an additional RF field Hj . The magnetization vector then 
10 comprises a rotating magnetic dipole m in the transverse plane. The NMR image may 
be reconstructed by sampling the external field from a series of RF dipoles. The discrete 
signals are recorded by a detector array over the sample space comprising the xz-plane 
and the positive y-axis of FIGURE 11. 

15 T l and T 2 NMR Images 

The NMR active nuclei including protons posses both angular momentum and a 
magnetic moment. When nuclei are placed in a static magnetic field H 0 , they precess 
about the field at a frequency proportional to the magnitude of H 0 . The bulk 

20 magnetization M of each voxel comprises the vector sum of the magnetic moments from 
all of the nuclei in each voxel. If the precessing nuclei are then subjected to an additional 
rotating (RF) field H 4 , which is synchronous with the precession, their magnetic 
moments and thus M will precess about H a and rotate away from the primary field H 0 
by an angle <j) Hi in a coordinate frame which rotates at the Larmor frequency [3]. The 

25 precession about H t continues as long as Hj exists. The final value of 0 Hi then depends 
on the strength of , which determines the precession rate, and the time for which it is 
turned on. The nuclei absorb energy as they change their orientation. This is known as 
nuclear magnetic resonance (NMR). The temperature of the nuclei or nuclear spin 
system rises during absorption of energy. When the Hj field is removed, the spin 

30 system cools down until it is thermal equilibrium with its environment. The exponential 
relaxation of the spin system temperature to that of the surrounding lattice is called spin- 
lattice relaxation and has a time constant T x where a time constant is defined as the time 
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it takes for 63% of the relaxation to occur. The NMR signal may also decay because the 
nuclei initially in phase following the Hj pulse get out of alignment with each other or 
dephase by local interactions with the magnetic fields of neighbor nuclei. The dephasing 
of the NMR signal is due to differing precession rates effected by the local interactions 
5 and is described by an exponential time constant T 2 also known as spin/spin relaxation. 
The main source of NMR image ( also called magnetic resonance images (MRI)) 
contrast is T x and T 2 which depend on tissue types. 

In an embodiment of the present invention, a T x image is produced by a applying 
at least one pulse sequence that inverts the magnetization and records the relaxation, a 

10 technique called inversion recovery. For example, the RF receivers are switched on to 
follow the decay following the nuclear excitation comprising a Hj pulse. M z (t) 9 the 
time dependent bulk magnetization in the direction of the primary field H 0 (coordinates 
of FIGURE 1 1) is examined at a time t Q after an inverting pulse by applying another Hj 
pulse equivalent to a rotation by 90° after waiting the time t Q following the initial 

15 inversion. The 90° pulse puts the z magnetization M z (* 0 ) into the transverse plane for 
observation. Changing the waiting time t 0 allows for observation of M z (t) at different 
times during relaxation. 

In an embodiment of the present invention, a T 2 image is produced by a applying 
at least one pulse sequence that flips the magnetization vector into the transverse plane 

20 and records the transverse relaxation by producing at least one spin-echo. The T 2 image 
depends on the NMR signal decaying because the nuclei initially in phase following the 
Hj pulse get out of alignment with each other or dephase by local interactions with the 
magnetic fields of neighbor nuclei. The dephasing of the NMR signal is due to differing 
precession rates effected by the local interactions and is described by an exponential time 

25 constant T 2 . Another unwanted source of dephasing is due to inhomogeneities in the 
primary field H 0 across the image space that causes an addition contribution to the 
precession rates of the magnetic nuclei. The signal ones observes after flipping the 
magnetization into the transverse plane that includes dephasing from an inhomogeneous 
H 0 field is still known as the free induction decay (FID), but the is called T* relaxation. 

30 The T 2 relaxation may be recovered from a T* FID since the H 0 inhomogeneity is 
constant and may be reversed. The dephasing due to the static inhomogeneity of H 0 
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may be canceled out by applying a 180° pulse at time t 0 along the y R -axis as shown 
in FIGURE 1 1 where the T* relaxation occurs in time t 0 . After an additional time t 0 , 
the total time elapsed after the 90° H x pulse is 2t Q = t E which is the spin-echo time. At 
this time the dephasing due to the static inhomogeneity of H 0 is exactly canceled out; so, 
5 the relaxation is strictly due to those processes that create T 2 relaxation which is 
recorded. After each 1 80° pulse, another spin echo is formed. The envelop of the 
maximum amplitude of the spin echoes is the T 2 relaxation. In an embodiment, a pulse 
sequence to give the data for a T 2 image known as the Carr-Purcell-Meiboom-Gill 
(CPMG) sequence [3] comprises applying a 90° pulse along the x R -axis followed by a 
10 series of 1 80° pulses along the y R -axis at times t Q + 2nt 0 where n is an integer including 
zero. 

In other words, when nuclei are placed in a static magnetic field H 0 and then 
subjected to an additional rotating (RF) field H 1 3 which is synchronous with their 
precession, M will precess about H a and rotate away from the primary field H 0 by an 
15 angle $ Hi . The magnitude of M is a maximum initially and decays with time. This 

occurs by emission of the same multipolarity radiation that it absorbed and by transfer of 
energy to the surrounding lattice. The intensity of the radiation is a function of M and 
the coordinate position relative to the RF emitting voxel, hi the present invention, the 
measurement of the intensity of the RF signal is performed over time and space 

20 following T x and/or T 2 encoding pulses. The signal as a function of time at a given 
detector position is Fourier transformed to give the components each having an 
amplitude and a unique phase. Each set of components of the NMR signal over the 
sample space due to a given voxel is determined from the phase data and the detector 
positions. The location of each voxel is determined through the spatial variations of the 

25 intensity of the transverse NMR field of the set of components associated by phase. 



4D-MRI System 



An NMR apparatus used to generate, and measure the secondary field and 
30 reconstruct the image is shown in FIGURE 8, and the corresponding coordinate system 
is shown in FIGURE 10. The apparatus comprises 1.) a magnet including a 
superconducting magnet to magnetize a volume of an object or tissue to be imaged, 2.) a 
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means (magnetometer) to determine the primary or magnetizing flux over the image 
space in the absence of the object to be imaged. In one embodiment, NMR of a proton 
containing homogeneous phantom such as water is determined on a point by point basis 
to map the primary field, or the magnetic flux at multiple points is obtained 
5 simultaneously using the reconstruction algorithm described herein, 3.) a radio frequency 
(RF) generator and transmitter including an antennae such as saddle coils to excite the 
protons of the magnetized volume, 4.) a means including an antennae coil to sample the 
dipole component (z-component in terms of Eq. (1)) of the RF secondary magnetic field 
at the Nyquist rate in time over the proton free induction decay, 5.) a detector array of 

10 elements of the means to sample the time signals including an antennae array which is 
selectively responsive to the dipole component (z-component of Eq. (1)) of the RF 
magnetic field of the magnetic moments of the protons which are aligned along the 
transverse axis to ideally point sample the secondary magnetic field at the Nyquist rate 
over the spatial dimensions which uniquely determine the NMR image which is 

15 reconstructed from the measurements, 6.) an analog to digital converter to digitize the RF 
signals, 7.) a time Fourier transform processor to convert the signal at each detector over 
time into its Fourier components, 8.) a processor to associate the Fourier components due 
to each voxel by phase, 9.) a reconstruction algorithm processor including a Fourier 
transform processor to convert each set of components into a voxel location of the bulk 

20 magnetization in the image space in parallel over all the voxels to form an NMR image, 
and 10.) image processors and a display such that the NMR image can be rotated in 
space to be displayed from any perspective as a three dimensional or two dimensional 
(tomographic) image. 

25 Magnetizing Field 

In an embodiment of the present invention, the applied magnetizing field which 
permeates the object to be imaged including tissue is confined only to that region which 
is to be imaged. The confined field limits the source of signal only to the volume of 
30 interest; thus, the volume to be reconstructed is limited to the magnetized volume which 
sets a limit to the computation required, and eliminates end effects of signal originating 
outside of the edges of the detector array. In the NMR case, the field having a steep 
gradient at the edges limits the imaged region by providing a range of Larmor 
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frequencies wherein the data is comprises a narrow frequency band at a desired Larmor 
frequency. 

Detector Array 

5 

An embodiment of the NMR imager of the present invention comprises a detector 
array of multiple detector elements which are arranged in a plane. The array may be a 
two dimensional detector array which is translated over the third dimension during the 
scan, or it may be a three dimensional detector array. The individual detectors of the 

10 array may respond to a single component of the secondary magnetic field which is 

produced by the magnetized object including tissue where the component of the field to 
which the detector is responsive determines the geometric system function which is used 
in the reconstruction algorithm discussed in the Reconstruction Algorithm Section. The 
detectors ideally point sample the secondary magnetic field at the Nyquist rate over the 

1 5 spatial dimensions which uniquely determine the NMR image which is reconstructed 
from the measurements. 

Small antennas may measure the RF signals as point samples without significant 
decrease in the signal to noise ratio relative to large antennas by using impedance 
matching while minimizing resistive losses by using superconducting reactance 

20 elements, for example. In an embodiment, cross talk between antennas is ameliorated or 
eliminated by time multiplexing the signal detection over the array of antennas. The RF 
field may be sampled synchronously or the field may be sampled at known times so that 
the phase at any given detector may be related to that at any other detector. 

Micromagnetic field sensors that are used to detect the primary field in the 

25 absence of the object to be scanned include NMR detectors and superconducting 

quantum interference devices (SQUIDS). Additional devices have been developed that 
are based on galvanometric effects due to the Lorentz force on charge carriers. In 
specific device configurations and operating conditions, the various galvanomagnetic 
effects (Hall voltage, Lorentz deflection, magnetoresistive, and magnetoconcentration) 

30 emerge. Semiconductor magnetic field sensors include MAGFETs, magnetotransitors 
(MT), Van der Pauw devices, integrated bulk Hall devices including the vertical MT 
(VMT), and the lateral MT (LMT), silicon on sapphire (SOS) and CMOS 
magnetodiodes, the magnetounijunction transistor (MUJT), and the carrier domain 
magnetometer (CDM), magnetic avalanche transistors (MAT), optoelectronic magnetic 
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field sensors, and magnetoresistive magnetic field sensors. In the case of NMR 
measurement of the secondary field (and/or the primary field), the detector array 
comprises RF antennas described in the NMR Primary Magnet, Gradient Magnets, RF 
Generator, RF Transmitter, and RF Receiver Section. 

5 

Scanning Methods 

The NMR scan performed on the object to be imaged including a human 
comprises the following steps: 

10 

• The magnetic moments of nuclei including protons of the object to be imaged that are 

aligned by the primary field are further aligned by a radio frequency (RF) pulse or 
series of pulses. 

15 • The strength and duration of the rotating H t (RF) field that is resonant with the protons 
of the magnetized volume and is oriented perpendicularly to the direction of the 
magnetizing field is applied such that the final precession angle of the magnetization 
is 90 ° (0 Hi - 90° ) such that the RF dipole is transverse to the primary magnetizing 
field and perpendicular to the RF magnetic field detector. 

20 

• NMR pulse sequences which provide the signals for a T x or T 2 image maybe applied. 

For example a 90 ° pulse may be followed by a series of 180 ° pulses. One sequence 
is the Carr-Purcell-Meiboom-Gill (CPMG) sequence [3]. 

25 • The free induction decay signals are recorded. 

• The time dependent signals are Fourier transformed to give the intensity and phase of 

each component. The NMR signal of each voxel at any given detector gives rise to a 
Fourier component with a unique phase angle relative to the Fourier component of 
30 any other voxel of the phantom at that detector. 

• The matrix of Fourier components that correspond to the NMR signal of a given voxel 

over the detectors is determined. This may be achieved by using a first component 
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having a phase angle and calculating the phase angle as a function of spatial position 
of the first detector relative to any other detector and identifying the component at 
each detector having the calculated phase angle. The matrices are determined for all 
of the voxels. The measurements of the spatial variations of the transverse RF field 
5 of a given matrix is used to determine the coordinate location of each voxel. Thus, 

each matrix of components associated by phase comprises the intensity variation over 
the sample space of the RF field of the bulk magnetization M of each voxel. 

• The Fourier transform algorithm given in the Fourier Transform Reconstruction 

10 Algorithm Section is performed on each set of components over the detector array to 
map each bulk magnetization M corresponding to a voxel to a spatial location over 
the image space. 

• In one embodiment, the point spread of the reconstructed voxel is corrected by 

1 5 assigning one voxel above a certain threshold with the bulk magnetization M . The 
other voxels are assigned a zero value. Ia this case, the reconstruction is digital 
versus analog. In terms of the signal to noise ratio, the superiority of digital over 
analog is generally known to those skilled in the art of signal processing. 

20 •In the limit with sufficient phase resolution, each volume element is reconstructed 

independently in parallel with all other volume elements such that the scan time is no 
greater than the nuclear free induction decay (FID) time. 

NMR Primary Magnet, Gradient Magnets, RF Generator, RF Transmitter, and RF 

25 

Receiver 

The primary magnet is that typical of a magnetic resonance imaging apparatus. 
The primary magnetizing field is in the z-direction as shown in FIGURE 10, usually 
30 coaxial with the patient's body. However, in an another embodiment the primary 
magnetizing field is in the direction which is perpendicular to the patient's body. An 
open primary magnet design is preferable such as a Helmholtz coil design to 
accommodate the three dimensional RF antennae array. A further field is applied to have 
a gradient in the z-direction (FIGURE 10). This provides a unique field in a chosen 
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volume of the patient including a cross-sectional slice of the patient. Li another 
embodiment, the NMR (Larmor) frequency from each voxel is determined by the 
magnitude of the magnetizing flux at the voxel position. The primary magnetic field has 
a sharp gradient at the edges of the image space wherein the Larmor frequencies outside 
5 of a selected range are rejected as arising from voxels outside of the image volume. 

A rotating RF field H t , of frequency chosen to cause resonance in the slice is 
then applied. Thus, only the molecules in the slice resonate. The resonance signal from 
the slice can then be detected. NMR pulse sequences which provide the signals for a T x 
or T 2 image may be applied. For example a 90 ° pulse may be followed by a series of 
10 180 ° pulses. One sequence is the Carr-Purcell-Meiboom-Gill (CPMG) sequence [3]. 

An embodiment of the NMR imager of the present invention comprises a RF 
generator 14 and 18 and RF transmitter coils 9 shown in FIGURE 8 which provide the 
rotating H t (RF) field that is resonant with the protons of the magnetized volume and is 
oriented perpendicularly to the direction of the magnetizing field as shown in FIGURE 
15 11. When the precessing nuclei are subjected to the additional rotating (RF) field H x , 
which is synchronous with the precession, their magnetic moments and thus M 
precesses about Hj and rotate away from the primary field H 0 by an angle 0 Hj in a 

coordinate frame which rotates at the Larmor frequency [3]. The precession about H 1 
continues as long as Hj exists. The final value of <j> Hi then depends on the strength of 

20 Hj , which determines the precession rate, and the time for which it is turned on. In an 
embodiment, the strength and duration of Hj is such that <j) Ui = 90° such that the dipole 
is oriented in the x R y R -plane of FIGURE 1 1 . A rotating coordinate frame is 
traditionally used to explain the physics of NMR [3]. Thus, is terms of the traditional 
NMR coordinate designation as described by Patz [3] and Hounsfield [7] shown in 

25 FIGURE 11, the x R -axis and y R -axis rotate about the primary field H 0 oriented parallel 
to the z R -axis at the Larmor frequency relative to the stationary NMR coordinate system 
(x, y, z) shown in FIGURE 10. Thus, the additional rotating (RF) field H t and the 
transverse RF magnetic field are stationary in the rotating NMR coordinate system, but 
rotate at the Larmor frequency in the stationary NMR coordinate system. Both the 

30 rotating and stationary coordinate systems are shown in FIGURE 11. In an embodiment, 
the rotating H t (RF) field is along the x R -axis and rotates the magnetization vector by 
the angle $ Hi = 90° to into the x R y R -plane. In terms of the reconstruction by the present 
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invention, the system function of the field corresponding to the rotated magnetization 
vector is equivalent to that of the detection of the z-component of a z-oriented dipole. 
The y-axis is the unique axis of the NMR system shown in FIGURE 8 and the stationary 
NMR coordinates shown in FIGURES 10 and 11. This axis corresponds to the unique 
5 axis of the Fourier Transform Reconstruction Algorithm Section, the z-axis shown in 
FIGURE 9, and the magnetization axis is the z-axis. The magnetization axis of the NMR 
rotating coordinates [3,7] is the z R -axis shown in FIGURE 1 1 . The measured transverse 
RF magnetic dipole oriented in the xy-plane is periodically parallel to the y-axis and 
rotates about the primary field H 0 oriented parallel to the z-axis at the Larmor frequency. 

10 The RF transmitter includes saddle coils. The RF receiver is a three dimensional 

array, or two arrays 180° from each other, or two orthogonal pairs of arrays, wherein 
each member of a pair is 180° from each other. The detector ideally point samples the 
RF field at the Nyquist rate described in the Reconstruction Algorithm Section. In one 
embodiment shown in FIGURE 1 1, each detector 20a is a coil antennae perpendicular to 

15 the y-axis and selective to the y-component of the RF secondary magnetic field arising 
from each voxel 14a magnetized in the z-axis direction. This two dimensional array is 
translated along the y-axis during a scan where readings of the secondary magnetic field 
are obtained as a function of the translation. In another embodiment shown in FIGURE 
12, the array 401 is three dimensional comprising multiple parallel two dimensional 

20 arrays 402 wherein each two dimensional array has a plurality of antennae coils 403 
which detect the RF field from each magnetized voxel 404. 

In another embodiment, the rotating Hj (RF) field shown in FIGURE 1 1 is along 
the y R -axis and rotates the magnetization vector by the angle $ Hi - 90° to project into 
the x R y R -plane. In terms of the reconstruction by the present invention, the system 

25 function of the field corresponding to the rotated magnetization vector is equivalent to 
that of the detection of the z-component of a z-oriented dipole. The x-axis is the unique 
axis of the NMR system shown in FIGURE 8 and the stationary NMR coordinates 
shown in FIGURES 10 and 11. This axis corresponds to the unique axis of the Fourier 
Transform Reconstruction Algorithm Section, the z-axis shown in FIGURE 9, and the 

30 magnetization axis is the z-axis. The magnetization axis of the NMR rotating 

coordinates [3,7] is the z R -axis shown in FIGURE 1 1 . The detector array receives 
similar signals as those produced by the transmitter, and both are perpendicular to the x- 
axis. Each detector is a coil antennae perpendicular to the x-axis and selective to the x- 
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component of the RF secondary magnetic field arising from voxels magnetized in the z- 
axis direction. This two dimensional array is translated along the x-axis during a scan 
where readings of the secondary magnetic field are obtained as a function of the 
translation. Li another embodiment, the array is three dimensional comprising multiple 
5 parallel two dimensional arrays 402 as shown in FIGURE 12 with the exception that the 
array 401 is perpendicular to the x-axis. The measured transverse RF magnetic dipole 
oriented in the xy-plane is periodically parallel to the x-axis and rotates about the 
primary field H 0 oriented parallel to the z-axis at the Larmor frequency. 

A linear combination of the cases of rotation of the RF magnetization along the 

10 x R -axis and the y R -axis is within the scope of the present invention and may be adopted 
in a manner straightforward to those skilled in the NMR art to apply the fields and 
detectors described for this invention. 

The method described herein may be performed on a suitable NMR examining 
apparatus such as that shown in simplified form in FIGURE 8, and the corresponding 

15 coordinate system is shown in FIGURE 10. Illustrated schematically only are coils 6, 
which provide B 0 , the steady primary field 7, which provide G x , the field gradient in the 
x-axis direction as shown, 8 which provide the G y , the field gradient in the y-axis 
direction as shown, 9 which provide the RF field, and 10 5 which provide G z , the field 
gradient in the z-axis direction as shown. The coils are driven by B 0 , G x , G y , RF, and 

20 G z drive amplifiers 11, 12, 13, 14, and 15 respectively, controlled by B 0 , G^, RF, and 

G z control circuits 16, 17, 18, and 19, respectively. These circuits can take suitable 
forms which will be well known to those with experience of NMR equipment and other 
apparatus using coil induced magnetic fields. The circuits are controlled by a central 
processing and control unit 20 to achieve the desired primary field, field gradients, and 

25 RF field to rotate the magnetization vector such that it is perpendicular to coil 9. 

The RF coils may be two saddle shaped coils 9 which are driven in parallel to 
provide the rotating RF field. The FID signals sensed are received in this example by the 
three dimensional array of RF coils 30 and are amplified by an RF amplifier 21 before 
being applied to signal handling circuits 22. The three dimensional detector array 30 is 

30 shown in more detail in FIGURE 12 as the three dimensional detector array 401. The 
circuits 22 are arranged to make any appropriate calibrations and corrections but 
essentially transmit the signals to the processing circuits to provide the required 
representation of the examined volume. These circuits can conveniently be combined 
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with the circuits which control the primary field, field gradients, and RF field and thus 
are included in the circuits indicated at 20. The picture thus obtained is viewed on a 
display 23, such as a television monitor, and this may include inputs and other 
peripherals 24 for the provision of commands and instructions to the machine, or other 
5 forms of output. The display is not limited and includes any medium of conveing the 
image. Examples of displays, but not limited to, include printers, cathode ray tube 
displays, liquid crystal displays, plasma screens, three dimensional modelers, 
holographic displays, laser monitors, and projection monitors. 

The apparatus also includes field measurement and error signal circuits 25 which 
10 receive signals via amplifiers 26 from field probes X x , X 2 , Y l , and Y 2 shown. 

The patient 27 is inserted in the tubular former of G x and G y coils 7, 8 and is 

supported there by a suitable couch or other supporting means. Such supports may be 
readily provided in any suitable form. 

The coils 7, 8 are two sets of coils axially displaced, each set comprising two 

15 pairs of saddle coils the pair 7 being at 90° to the pair 8. These coils are themselves 
inserted into the central aperture in B 0 coils 6 which in an embodiment are wound in 
four parts connected in series to provide an approximately circular configuration which 
is well known to be desirable for production of a uniform field. Further details of the 
coil winding will not be given since suitable coils can readily be devised, by those with 

20 the appropriate skills, to provide the fields required. 

The appropriate stores provide the amplitude and duration signals which are 
converted to analog form in digital to analog converters (DAC ? s) and applied to 
respective coil drive circuits x, y, z, RF. The respective drive circuits, which can take 
any form well known for driving field coils, provide the specified currents to the 

25 appropriate coil for the specified duration. The apparatus and circuits described so far 
may be adopted to provide different gradients and RF fields, by appropriately adjusting 
the stored sequences and profile data. Similarly other known NMR apparatus which are 
capable of applying a steady magnetic field, a pulsed RF field, and G x , G y , and G z field 

gradients to a body, may be adopted in a manner straightforward to those skilled in the 
30 NMR art to apply the fields described for this invention. 

Gradients may be applied in any direction to further enhance the image 
reconstruction or image quality by methods known to those skilled in the art. In an 
embodiment shown in FIGURE 8, the system components which provide imaging 
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enhancing gradients are 7, which provide G x , the field gradient in the x-axis direction as 
shown, 8 which provide the G y , the field gradient in the y-axis direction as shown, and 
10, which provide G z , the field gradient in the z-axis direction as shown. For example, a 
magnetic field gradient along an axis in the xy-plane of the z-directed primary field H 0 
5 is applied to produce a phase dispersion as a function of the distance along the axis. In 
an embodiment, the phase due to the gradient may be linear and may be larger than the 
phase of any Fourier component in the absence of a gradient. During reconstruction, the 
phase information due to the gradient may be applied to refine the assignment of the 
position of each voxel with respect to the distance along the gradient axis by methods 

1 0 known to those skilled in the art. 

With the basic signal handling system of the present invention, the FID signals 
from the signal sensing coils of the detector array 30 (shown in more detail in FIGURE 
12 as coils 403 of the three dimensional detector array 401) are amplified in an RF 
amplifier and applied via an analog to digital converter (ADC) to a store such as a 

1 5 random access memory (RAM). The data is then processed according to the procedure 
given in the Reconstruction Algorithm Section. 



Phase Angle and Associated Set of Fourier Components 

20 

In an embodiment, the NMR signal at each detector as a function of time is 
processed by a method such as a Fourier transform operation to give a plurality of 
Fourier components each having an intensity, a phase angle, and the same frequency. 
The NMR signal of each voxel at any given detector gives rise to a Fourier component 

25 with a unique phase angle relative to the Fourier component of any other voxel of the 
phantom at that detector. The set of Fourier components that correspond to the NMR 
signal of a given voxel over the detectors is determined. This may be achieved by using 
a component with a defined standard phase relative to a first detector. Over the array of 
sampled points, the phase angle of components at each detector may be converted to the 

30 corresponding phase at the position of the component having the standard phase. The 
components that have the standard phase are associated into a set of components 
comprising the components from a given voxel over the sample space. 
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Li an embodiment, the phases of the components from the voxel are determined 
at the other detectors using the relative position of other detectors of the sample space 
relative to the first detector and the standard phase of the component at the first detector. 
In other words, the set of associated components may also be determined by 1.) 
5 identifying a first component having a characterizable phase angle, 2.) calculating the 
phase angle as a function of spatial position of the other detectors relative to the first 
using the phase angle of the component at the first detector, and 3.) identifying the 
component at each detector having the calculated phase angle. These components 
identified from the data are associated to form the set. 

10 In the present invention, the phase variation of the NMR signal is provided by the 

combination of 1 .) the angle 9 suspended between the direction of the detector and the 
radial vector, the vector from the dipole to the detector, and 2.) the angle $ due to a 
separation distance r between a voxel and a detector given by the wavenumber of the RF 
field k times r . The distance r of the detector from the voxel gives rise to the phase 

15 term e~ lh ' of the component of the detected RF signal where k is the wavenumber of the 
NMR signal. The harmonic oscillation of each RF dipole is equivalent to the dipole 
rotating in the transverse plane. The detector is responsive to a component in this plane. 
At a point in time, each RF dipole is directed at an angle 0 relative to the direction of 
detection of the detector. The angle 9 of the RF dipole relative to the direction of 

20 detection axis of the detector gives rise to a phase angle term e~ lQ . In a preferred 

embodiment, the sum of the phase angles, kr and 9 , are unique for each voxel at each 
detector. The position of each detector relative to a different detector may be used to 
calculated the phase angle of the second relative to the first. This may be repeated over 
all of the detectors to give the set of intensities of the NMR signal over the sample space 

25 due to each voxel. 

In an embodiment, the voxels that are on the same axis with a given detector or a 
plurality of detectors that align vertically with one or more voxels is determined. The 
phase angle is then only a function of the distance r of the detector from the voxel. Such 
components may be identified by the presence of at least one component of a phase 

30 given by 0 + kAr where <j) x is the phase angle of the first component and Ar is the 
distance of the second detector relative to the first of the detector. This detector(s) is 
defined as having a zero angle 9 . In this case, the phase angle 0 2 of a second component 
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at a second detector aligned on the detection axis corresponding to a first component at a 
first detector having a phase angle $ x is given by 



From the components identified as coaxial with a detector, the phases of the components 
from the voxel are determined at the other detectors using the relative position of other 
detectors of the sample space relative to the first detector and <^ of the component at the 
first detector. 

FIGURE 14 shows the stationary coordinate system (x, y, z) of the NMR detector 
corresponding to the coordinate system of FIGURE 9 and FIGURE 10 of a primary field 
H 0 oriented parallel to the z R -axis and the z-axis and the corresponding transverse RF 
magnetic dipole oriented in the x R y R -plane and periodically parallel to the y-axis 
wherein the spatial variation of the RF y-component of the flux due to the RF dipole is 
measured at a detector according to Eq. (1). FIGURE 14 further shows the distances and 
angles between a voxel 801 linear to a first detector 802, a second nonlinear voxel 800, 
and a second detector 803 in accordance with the invention. The NMR signal of each 
voxel at any given detector gives rise to a Fourier component with a unique phase angle 
relative to the Fourier component of any other voxel of the phantom at that detector. In 
the present invention, the unique phase variation of the NMR signal is provided by the 
combination of 1.) the angle 0 suspended between the direction of the detector and the 
radial vector, the vector from the dipole to the detector, and 2.) the angle <p due to a 
separation distance r between a voxel and a detector. The phase <j) due to a separation 
distance between a voxel and a detector of r is given by the wavenumber of the RF field 
k times r . 



where the wavenumber k is given by 
k = — (20) 

In terms of the NMR coordinates of the detector shown in FIGURE 1 1 and FIGURE 14, 
the phase angel <j) is given by 



02 = 0i + 



(18) 



(19) 




(21) 
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where Ax, Ay, and Az is the separation distance between a voxel and a detector in the 
x , y , and z directions, respectively. In the case that a two dimensional slice is imaged 
in the transverse plane, the phase angel (j) is given by 

<j>=kr = ^(Ax) 2 + (Ay) 2 (22) 

5 All voxels which have the same radial distance from a given detector will have 

the same phase given by Eq. (19). However, at each case that the radial distances r of 
two voxels to a given detector are equivalent, a unique angle 9 suspended between the 
direction of the detector and the vector along r exists. In an embodiment, the RF dipoles 
of all voxels are time synchronous. The NMR pulsing and detection may be 

10 synchronized. At least one standard dipole may be used as reference phase to set the 
phases of the RF dipoles of the phantom. All voxels that are not collinear with the axis 
of a detector will posses an angle 9 suspended between the direction of the detector and 
the radial vector, the vector from the dipole to the detector. Since the RF dipoles are 
time synchronous, at each point in time, this angle corresponds to a phase angle 9 of the 

15 RF signal from each voxel at each detector. In the case that a two dimensional slice is 
imaged in the transverse plane, the phase angle 9 is given by 

9 = sm — = sm , 0 = sm . (23) 

r tI(ax) + (Ay) V ( kAx ) + ( kA y) 

The total phase angle <j) T of the RF signal from each voxel at each detector is given as the 
sum of <j> and 9 . 

20 <j> T =(j) + 6 (24) 

The uniqueness of the phase angle a Fourier component of the NMR signal of 
each voxel at any given detector relative to a Fourier component of any other voxel of 
the phantom at that detector is demonstrated by the following cases according to the 
coordinates shown in FIGURE 14: 

25 

Case 1 Ax = 0.1 cm Ay = 0.1 cm 

0=sin - i - f — 2i — = 4>=^V(o.i) 2 +(o.i) 2 

V(o.i) 2 +(0.1) 2 X 
9 = 0.78 rad (j> = 5.9 X10 -4 rad 

30 <\> T =0 + 0=0.78 + 5.9X10^ rad= 0.7806 rad 
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Case 2 Ax = 0.1 cm Ay = 20 cm 

e= sin -' =L_ <^ = — V(0.l) 2 + (20) 2 

V(0.1) 2 +(20) 2 ^ ^ ^ V ^ 

5 0=5XlO~ 3 rad <j> = 8.37 X10" 2 rad 

0 r =0 + 0=5 X10~ 3 + 8.377 X10" 2 rarf = 8.88X10" 2 rarf 

Case 3 Ax = 20 cm Ay = 0.1 cm 

10 e^sin- 1 , , 2 ° , — V(20) 2 +(0.1) 2 

V(20) 2 +(0.l) 2 ^ A W 

0 = 1.56 rad <j> = 8.37 X10" 2 rad 

(j> T =9 + <j)=1.56 + 8377X10~ 2 rad = 1.643 rad 



15 



Case 4 Ax = 20 cm Ay = 20 cm 

e=shrl J 0 ,,, ^=^V(20) 2 + (20) 2 

V(20) 2 + (20) 2 & 

6 =0.785 rad =0.1184 rarf 

0 r =9 + <f> =0.785 + 0.1184raJ = 0.9034 raJ 

The condition for uniqueness of the phase 0 r of each voxel at a given detector 
20 depends on the unique sum given by Eq. (24). Since all voxels which have the same 

radial distance from a given detector will have the same phase (j> given by Eq. (19), but a 
unique angle 9 , each voxel will have a unique phase <j> T when the following condition is 
satisfied: 

6^ , the maximum angle 9 , is greater than 0^ , the maximum angle (j) , for any 

25 r. 

0,«>4L« (25) 
9^ corresponds to x^ and y^ , the minimum separation of the voxel and the detector 
in the x and y directions, respectively, and is given by Eq. (23). 
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6^ = sin' 1 , *™ 2 (26) 

$ max corresponds to x min , the minimum separation of the voxel and the detector in the x 
direction and >> max ? the maximum separation of the voxel and the detector in the y 
direction, and is given by Eq. (22). 

5 =y ^(^y+iy^y (27) 

Substitution of Eq. (26) and Eq. (27) into Eq. (25) gives the condition for uniqueness of 
the phase angle <p T . 

1° For » ^min . Eq. (29) gives 

\( *\ V >f y - (30) 

For = ^ , Eq.(30) gives 

X>27tj2y^ (31) 
^>8.9 ymax (32) 
15 In terms of frequency /, Eq. (32) gives 

f< TZ~ (33) 

where c is the speed of light. For y max = 20 cm, the condition for uniqueness of the 
phase angle <j) T is / =< 168 MHz. 

Consider a component of a voxel at a first detector with a phase <^ defined as a 

20 standard phase. The corresponding phases of the components from the same voxel are 
determined at other detectors using the relative position of the other detectors of the 
sample space relative to first detector. In the case that a voxel is collinear with the first 
detector, the phase angle 0=0. The phase angle <j) 2 of the corresponding component at 
a second collinear detector can be calculated from the phase angle <^ of a first using the 

25 detector separation along the collinear axis. In the case that a two dimensional slice is 



WO 02/16956 PCT/US01/25954 

41 

imaged and the second detector is displaced Ay in the y direction of the sample space 
relative to the first detector with a phase angle <£ 1? the angle <f> 2 is given by 

<j> 2 = V(W) 2 + (kAyf = V(W f + (34) 
Consider the phase angle § T2 of a corresponding component at a second detector 
5 that is nonlinear with the first detector. The phase angle 0 r2 of the second component 
due to the same voxel may be calculated from § v In this case a contribution exists from 
1 .) the angle 6 suspended between the direction of the detector and the radial vector, the 
vector from the dipole to the detector, and 2.) the angle $> 2 due to a separation distance r 
between a voxel and a detector given by the wavenumber of the RF field k times r . In 
10 the case that a two dimensional slice is imaged and the second detector is displaced Ax 1 
and Ay in the x and y directions of the sample space, respectively, relative to the first 
detector with a phase angle <f> l? the phase 9 is given by 



« . _i Ax' . _i Ax' . _i kAx } 

6 = sin 1 = sm 1 j = sin 1 i 

* V( Ajc ') 2 + (Av f ) 2 s( kAx ' f + ( kA yJ 

. _! kAx" 



(35) 



= sm 



V(^Ax') 2 + (My) 2 + ^ 2 

The angle <j> 2 is given by 
15 <j> 2 = ) 2 + (My) 2 + (kAy) 2 = ^(kAx'f + (My ) 2 + 0 2 (36) 

The phase angle </> r2 of the second component due to the same voxel is given by the sum 
ofEq. (35)andEq. (36). 

fe^e+^sfaT 1 , +V(^) 2 + (My) 2 +0 1 2 (37) 

V(MxO +(My)P + ^ 2 

The component at the second detector which has phase angle <j) T2 is associated with the 
20 component from which <f> T2 was calculated. All phase angles <f> T2 at all other detectors 
are determined from each standard phase angle $ r This is repeated for all using those 
detectors which are determined to are collinear with a voxel to give the sets of 
components each comprising the spatial variation of the RF field over the sample space 
due to each voxel. 
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Reconstruction Algorithm 
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When a static magnetic field H 0 with lines in the direction of the z-axis is applied 
to an object comprising a material containing nuclei such as protons that possess 
5 magnetic moments, the field magnetizes the material. As a result a secondary field 
superposes the applied field as shown in FIGURE 9. In the applied magnetic field, the 
magnetic moments of each nuclei precesses about the applied magnetic field. However, 
the magnetization of any one nucleus is not observed from the macroscopic sample. 
Rather the vector sum of the dipole moments from all magnetic nuclei in the sample is 

10 observed. This bulk magnetization is denoted by the vector M . Li thermal equilibrium 
with the primary field H 0 , the bulk magnetization M is parallel to H 0 . The volume to 
be imaged is divided into volume elements called voxels, and the magnetized voxel 302 
shown in FIGURE 9 with bulk magnetization M is modeled as a magnetic dipole m . 
Consider the case wherein data comprising the z-component of the magnetic field of a 

15 dipole oriented in the z-direction is acquired by detectors 301 in the three dimensional 
sample space comprising the xy-plane and the positive z-axis as shown in FIGURE 9. 
The magnetic moment, m z , of each voxel within the phantom is a magnetic dipole. And, 
the phantom can be considered to be a three-dimensional array of magnetic dipoles. At 
any point extrinsic to the phantom, the z-component of the secondary flux, B , from any 

20 single voxel is 



2z 2 ~x l ~y 2 

B =^7T-T— (38) 
(x +y +z ) 

where x, y, and z are the distances from the center of the voxel to the sampling point. It 
is shown in APPENDICES I-IV that no geometric distribution of magnetic dipoles can 
give rise to Eq. (1). Therefore, the flux of each magnetic dipole (voxel contribution) 

25 serves as a basis element to form a unique reconstruction of the array of dipoles which 
comprise the bulk magnetization map or NMR image of the phantom. 

Eq. (1) is a system function which gives the magnetic flux output in response to a 
magnetic dipole input at the origin. The phantom is an array of spatially advanced and 
delayed dipoles weighted according to the bulk magnetization of each voxel; this is the 

30 input function. The secondary flux is the superposition of spatially advanced and 
delayed flux, according to Eq. (1); this is the output function. Thus, the response of 
space to a magnetized phantom is given by the convolution of Eq. (1) with the series of 
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weighted, spatially advanced and delayed dipoles representing the bulk magnetization 
map of the phantom. The discrete signals are recorded by a detector array over the 
sample space comprising the xy-plane and the positive z-axis. 

hi an embodiment of the present invention, the magnetization vector is rotated 
5 into the transverse plane by an additional RF field H x . The magnetization vector then 
comprises a rotating magnetic dipole m in the transverse plane. The NMR image may 
be reconstructed by sampling the external field from a series of RF dipoles rather than 
that from a series of static dipoles. In this case, the geometric system function is given 
also given by Eq. (1), the function of the z-component of the flux B from a z-oriented 
10 dipole at any point extrinsic to the image space, from any single voxel. The geometric 
system function corresponding to a dipole which determines the spatial intensity 
• variations of the RF field is a band-pass for k p = k z as shown in APPENDIX V. In the 

limit, each volume element is reconstructed independently in parallel with all other 
volume elements such that the scan time is no greater than the nuclear free induction 
15 decay (FID) time. 

In an embodiment, the strength and duration of the rotating H t (RF) field that is 
resonant with the protons of the magnetized volume and is oriented perpendicularly to 
the direction of the magnetizing field is applied such that the final precession angle of the 
magnetization is 90° (0 Hi = 90° ) such that the RF dipole is transverse to the static 

20 magnetizing field and perpendicular to the RF magnetic field detector 20a as shown in 
FIGURE 1 1 . The rotating NMR coordinates and the stationary NMR coordinates are 
both shown in FIGURE 11. According to Eq. (1), the signal as a function of time which 
is Fourier transformed arises from each transverse RF magnetic dipole oriented in the 
x R y R -plane (rotating NMR coordinates) which is periodically parallel to the y-axis and 

25 rotates about the primary field H 0 oriented parallel to the z-axis at the Larmor frequency 
(stationary NMR coordinates). For the stationary NMR case, the y-axis corresponds to 
the z-axis of Eq. (1), and m z , the magnetic moment along the z-axis, of Eq. (1) 
corresponds to the bulk magnetization M of each voxel. 

In addition, each Fourier component comprises an additional part that 

30 corresponds to an RF magnetic dipole which is periodically parallel to the x-axis. This 
part also rotates about the primary field H 0 oriented parallel to the z-axis at the Larmor 
frequency wherein the x-axis of the stationary NMR coordinates corresponds to the x- 
axis of Eq. (1). The parts rotate at the same frequency but are orthogonal. The 
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orthogonality corresponds to a phase angle of the Fourier transform of 90°; thus, each 
Fourier component contains a real and an imaginary part. The magnitude of the 
amplitude of the signal due to the dipole oriented parallel to the y-axis (z-axis of Eq. (1)) 
exceeds that of the dipole oriented parallel to the x-axis (x-axis of Eq. (1)) which 
5 identifies this part of a Fourier component. In an embodiment, the spatial variation over 
the detector array of the part of the Fourier component with the maximum amplitude (the 
signal due to the dipole oriented parallel to the y-axis) is used to determine the 
coordinate location of each voxel using the geometric system function of the detection of 
the z-component of a z-oriented dipole (Eq. (1)) as given in the Reconstruction 

10 Algorithm Section. 

In another embodiment, the spatial variation over the detector array of the lesser 
magnitude orthogonal part of each Fourier component is used to determine the 
coordinate location of each voxel using the corresponding geometric system function as 
described in the Reconstruction Algorithm Section. The system function for the case of 

15 the minor orthogonal part in terms of the coordinates given in FIGURE 9, corresponds to 
the detection of the x-component of a z-oriented dipole. The geometric system function 
is given by the function of the x-component of the flux B from a z-oriented dipole at any 
point extrinsic to the image space, from any single voxel: 
, 3xz 

B = ™z ,2^2^ 2x5/2 (39) 

20 where x, y ? and z are the distances from the center of the voxel to the sampling point and 
m z is the magnetic moment along the z-axis corresponding to the bulk magnetization M 
of each voxel. 

In the reconstruction process described herein, the secondary field may be in the 
same or in a transverse orientation relative to the primary field. In each case, the 

25 orientation of the secondary dipole field and the measured secondary dipole field 

component are according to Eq. (1). The relationship of the coordinate systems of the 
present invention are shown in FIGURES 9, 10, and 11. The primary and secondary 
fields are parallel and stationary in FIGURE 9 versus transverse in the rotating and 
stationary NMR coordinates shown in FIGURES 10, 1 1, and 12. The relationship of the 

30 coordinate system of the reconstruction method of the NMR system used herein to that 
of the coordinate system of Eq. (1) shown in FIGURE 9 is according to FIGURE 1 1 . 
The magnetization axis of Eq (1) is the z-axis; whereas, the magnetization axis of the 
stationary NMR coordinates is the z-axis shown in FIGURE 10, and the magnetization 
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axis of the NMR rotating coordinates [3,7] is the z R -axis shown in FIGURE 1 1 . 
Regarding Eq. (1) and the reconstruction algorithm in the NMR case, the z-component of 
Eq. (1) is substituted with the stationary NMR y-component, the y-component is 
substituted with the stationary NMR z-component, and the x-component is substituted 
5 with the stationary NMR x-component. 

The reconstruction algorithm can be a reiterative, a matrix inversion, or a Fourier 
Transform algorithm. For all reconstruction algorithms, the volume to be imaged is 
divided into volume elements called voxels, and the magnetized voxel with magnetic 
moment m z is modeled as a magnetic dipole. In an embodiment, the matrix of Fourier 

1 0 components that correspond to the NMR signal of a given voxel over the detectors is 

determined. The matrices are determined for all of the voxels. The measurements of the 
spatial variations of the transverse RF field of a given matrix is used to determine the 
coordinate location of each voxel. Thus, each matrix of components associated by phase 
comprises the intensity variation over the sample space of the RF field of the bulk 

15 magnetization M of each voxel. The matrices are the input for the reconstruction 

algorithm. This procedure is performed either in parallel or series for each matrix. The 
bulk magnetization map (NMR image) is the superposition of the independent images, 
each of which corresponds to a given voxel. The superposition of images is plotted and 
displayed. 

20 An embodiment of a matrix inversion reconstruction algorithm comprises the 

steps of 1.) using the geometric system function to determine the spatial intensity 
variation of the transverse RF field over the detector array, 2.) inverting the 
corresponding matrix, and 3.) multiplying the signal over the detector array by the 
inverted matrix to give the voxel sources. For example, a matrix inversion 

25 reconstruction algorithm is to determine a coefficient for each voxel mathematically (Eq. 
(41)) or by calibration which when multiplied by the bulk magnetization M of each 
voxel is that voxel's contribution to the signal at a given detector with the corresponding 
unique phase at each detector. This is repeated for every detector, and those coefficients 
are used to determine a matrix which, when multiplied by a column vector of the bulk 

30 magnetization M values of the voxels, gives the voltage signals at the detectors. This 
matrix is inverted and stored in memory. Voltages as a function of time are recorded 
over the detector array. The signal as a function of time is Fourier transformed to give 
the Fourier components each having an amplitude and a unique phase. In an 
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embodiment, of a matrix inversion reconstruction algorithm, the components are 
multiplied by the inverse matrix, to generate the bulk magnetization M map. 

In another embodiment of a matrix inversion reconstruction algorithm, the matrix 
of Fourier components that correspond to the NMR signal of a given voxel over the 
5 detectors is determined. The matrices are determined for all of the voxels. The 
measurements of the spatial variations of the transverse RF field of a given matrix is 
used to determine the coordinate location of each voxel. Thus, each matrix of 
components associated by phase comprises the intensity variation over the sample space 
of the RF field of the bulk magnetization M of each voxel. Each matrix of components 

10 are multiplied by the inverse matrix, to generate the bulk magnetization M map. In one 
embodiment, the point spread of the reconstructed voxel is corrected by assigning one 
voxel above a certain threshold with the bulk magnetization M . The other voxels are 
assigned a zero value. This procedure may be repeated for all voxels. In the limit with 
sufficient phase resolution, each volume element is reconstructed independently in 

15 parallel with all other volume elements such that the scan time is no greater than the 
nuclear free induction decay (FED) time. The total bulk magnetization map (NMR 
image) is the superposition of the separate maps for each magnetic moment which is 
plotted and displayed. 

In an embodiment of a reiterative reconstruction algorithm, the geometric system 

20 functiomis used to determine the system of linear equations which gives the intensity, 
spatial variation, and phase of the RF field over the sample space. The signal as a 
function of time is Fourier transformed to give the Fourier components each having an 
amplitude and a unique phase. The system of linear equations gives the voltage and 
phase from each voxel at each sensor based on the bulk magnetization M value of each 

25 voxel and the position of the voxel relative to the sensor. In an embodiment, weighting 
coefficients are determined based on these equations. The coefficients may be 
determined mathematically (Eq. (41)). alternatively, they may be determined by 
calibration. The bulk magnetization M for each voxel is estimated, and the signals at 
each detector are calculated at each phase. In an embodiment, the bulk magnetization M 

30 value of each voxel times its weighting coefficient and its calculated phase at a given 
detector is compared to the measured voltage and phase. A correction is made to M of 
each voxel. This gives rise to a second, or recomputed, estimate for M of each voxel. 
The signal value from this second estimate is computed and corrections are made as 
previously described. This is repeated until the correction for each reiteration 
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approaches a predefined limit which serves to indicate that the reconstruction is within 
reasonable limits of error. This procedure is repeated for all sensors. The final bulk 
magnetization map is plotted and displayed. 

The general process of reconstruction by reiteration is shown according to the 
5 steps of FIGURE 13 (and is implemented in processor 20 in FIGURE 8). The image 
displayed according to the process 200 is merely a mapping of the bulk magnetization 
M of voxel sections of the object examined. The signal as a function of time is Fourier 
transformed to give the Fourier components each having an amplitude and a unique 
phase. The bulk magnetization M at a given location determines the intensity and 

10 spatial variation of the RF field over the sample space. Accordingly, signals produced 
by the RF magnetic sensors 1 10, in terms of volts as a function of phase, are a direct 
result of the bulk magnetization M of the voxel elements. Therefore, a reference voltage 
is generated at 210 from which the actual or measured sensor voltage at a given phase is 
subtracted at 220. The reference voltages are modeled by assuming a voltage with a 

15 unique phase from each voxel at each sensor. Therefore, the voltage and phase from 
each voxel is determined at a sensor according to a weighting based on the position of 
the voxel relative to the sensor. The weighting may be given by Eq. (1) in terms of the 
coordinates x, y, z of FIGURE 9. The phase may be determined as given in the Phase 
Angle and Associated Set of Fourier Components Section. The resulting modeled or 

20 calculated voltage signals are compared at step 220, providing a difference or A(x, y, z) 
signal, weighted at step 230 to produce a weighted difference signal, which is then added 
to the previously estimated bulk magnetization M value for each voxel element at step 
240. The resulting level, available in three dimensions corresponding to the axes x, y, 
and z, is selectively displayed on the display at step 250. Having adjusted the estimated 

25 bulk magnetization M for each voxel, the calculated bulk magnetization M is 

recalculated at step 260, the resulting estimated sensor voltage at each phase is then 
compared to the actual sensor voltage at each phase at step 220, the process 200 being 
repeated until the difference is reduced to a tolerable value. This procedure is repeated 
for all sensors. The final bulk magnetization map is plotted and displayed. 

30 In terms of the coordinates of Eq. (1) and FIGURE 9, the reconstruction 

algorithm using Fourier Transforms (FFT) involves exploiting the FFT to solve Eq. (41) 
given below. For the case that follows, data is acquired in the x, y, and z-directions, but 
in general, data is acquired over the dimensions which uniquely determine the bulk 
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magnetization map (NMR image). Also, the present analysis is for measuring the z- 
component of the magnetic field of a dipole oriented in the z-direction; however, the 
analysis applies to the other two orthogonal components where the geometric system 
function for the z-component of the z-oriented dipole is replaced by the geometric 
5 system function for the x or y-components of the magnetic field produced by the dipole 
where the geometric system function is defined below as the impulse response of the 
detector to the given component of the field of a dipole of given orientation. The sample 
space, or space over which the secondary RF field is measured, is defined in the present 
example as the three-dimensional space comprising the entire xy-plane and the positive 

10 z-axis, as shown in FIGURE 9. Other sample spaces are valid and each requires special 
consideration during the reconstruction as described below. The discrete voltages 
recorded from an infinite detector array in the xy-plane which is translated to infinity 
along the z-axis starting from the origin where the detector array is responsive to the z- 
component of the secondary magnetic field is given by Eq. (1), where the voltage at any 

15 point in space produced by dipoles advanced in the z-direction and advanced or delayed 
in the x and y-directions relative to the origin is given by the following Eq. (40), where 
the voltage response is C 0 times the secondary magnetic flux strength in the case shown 
in FIGURE 9. 

In a preferred embodiment, the data comprises the RF field in the transverse 
20 plane over the sample space shown in FIGURE 11. The relationship of the coordinate 
designations of Eq. (1) and FIGURE 9 versus the NMR system and the relationship 
between the magnetic moment of a static secondary field and the bulk magnetization M 
of the RF field is given in the Reconstruction Algorithm Section and in the Fourier 
Transform Reconstruction Algorithm which follows. 

25 

Fourier Transform Reconstruction Algorithm 

In terms of the coordinates of Eq. (1) and FIGURE 9, the volume to be imaged is 
30 divided into volume elements called voxels and the magnetized voxel is modeled as a 
magnetic dipole m t . Data comprising the z-component of the magnetic field of a dipole 
oriented in the z-direction is acquired in the three dimensional sample space comprising 
the xy-plane and the positive z-axis. 
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The phantom can be considered to be a three-dimensional array of magnetic 
dipoles. At any point extrinsic to the phantom, the z-component of the secondary flux, 
B , from any single voxel is 
2z 2 -x 2 -v 2 

B =m , 2 2 ; 5/2 (40) 
(x +y + z ) 

5 where x, y, and z are the distances from the center of the voxel to the sampling point. It 
is shown in APPENDICES I-IV that no geometric distribution of magnetic dipoles can 
give rise to Eq. (1). Therefore, the flux of each magnetic dipole (voxel contribution) 
forms a basis set for the flux of the array of dipoles which comprise the bulk 
magnetization map of the phantom. 

10 Eq. (1) is a system function which gives the magnetic flux output in response to a 

magnetic dipole input at the origin. The phantom is an array of spatially advanced and 
delayed dipoles weighted according to the bulk magnetization of each voxel; this is the 
input function. The secondary flux is the superposition of spatially advanced and 
delayed flux, according to Eq. (1); this is the output function. Thus, the response of 

15 space to a magnetized phantom is given by the convolution of Eq. (1) with the series of 
weighted, spatially advanced and delayed dipoles representing the bulk magnetization 
map of the phantom. The discrete voltages recorded by a detector array over the sample 
space comprising the xy-plane and the positive z-axis are 

«3=o n 2 =-i 2 iik ni =~ij2k H x _ n ik J + \y - J + [z + n 3 k 3 \ 

20 (41) 
where the flux magnetizing each voxel is unity, the volume element is made unity, and 
the voltage response of each detector is C 0 times the secondary magnetic flux strength. 
The variables of Eq. (40) are defined as follows: 

25 x m ,« 2 , m ^ magnetization of the voxel located at the 

position defined by the Dirac delta function, 
5(x —n x k^ y— njc 29 z - /7 3 /c 3 ) 



30 



k\ 9 k 2 , 



dipole spacing in the x, y, and z-directions, 
respectively 
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the dimensions of the image space, the total distance in the x, y, and z-directions 5 
respectively, for which the bulk magnetization of the voxels is nonzero 

the detector spacing in the x, y, and z-directions, respectively 

5 

The voltage signals recorded at the detector array over the sample space is 

oo oo co 

F[m 1 ,m 2 ,/n 3 ]= J) C o S [ x - m^y-m^z -m 3 s 3 \ 

/«3 = 0 ni2 =-oo m\— —oo 

+ ^2A \2 7 ?-x 2 -y 2 ~\ _ J 

Zj Z, lu T~2 1 7^ <S> %»un2,« 3 5 L X - n i k l>y- n 2 k 2> Z - 1 h k 3]\ 

»,=o )<2 =-/ 2 /2* «,=-V2tp: +y +z J J 

(42) 

Eq. (41) can be represented symbolically as follows 
10 s=C 0 [gx[h®f]xu(z)] (43) 

where C 0 is the proportionality constant between the signal voltage and the output flux 
strength, s is the discrete function of voltage signals recorded of the secondary flux over 
the sample space, g is the secondary magnetic flux sampling function given as follows: 

CO oo CO 

g = Yj 2 Yj S[x-m 1 s 1 ,y-m 2 s 2 ,z-}n 3 s 3 ] (44) 

m 3 =0 m 2 =- oo m l — -oo 

15 In Eq. (42), h is the system function which is also defined as the geometric system 
function given as follows: 
2z 2 -x 2 ~y 2 



\%2 



(45) 



(x 2 +y 2 +z 2 J 

The system function is the impulse response of the detector array where the external 
magnetizing flux is set equal to unity. In Eq. (42), / is the bulk magnetization function 
20 or NMR image function given as 

+ * 3 M + lj2k +lj2k 

S Z Xnu^]x-n l k i9 y-n 2 k 29 z-nJ^\ (46) 

7*3 = 0 ti 2 =-/ 2 / 2k n l =-l l /2k 

In Eq. (42), u(z) is the unitary z function which is one for positive z and zero otherwise. 
The function g discretizes the continuous voltage function V given by Eq. (40) which is 
the function h convolved with the function / then multiplied by the function u(z). The 
25 discrete voltages recorded over the sample space are used in a Fourier transform 
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algorithm derived herein to reconstruct the NMR image. Consider the function $ ofEq. 
(41) which is given as follows: 

s = (h®f)xu(z) (47) 
Eq. (46) is equivalent to the function h convolved with the function / then multiplied 
5 by the function u(z). The function S which is the Fourier transform of s is given as 
follows: 

S = (HxF)®U(k 2 ) (48) 
The function S is equivalent to the resultant function of the function H, the Fourier 
transform of the system function--/; , multiplied by the function F, the Fourier transform 
10 of the bulk magnetization function—/, convolved with the function U(lc z \ the Fourier 
transform of the unitary z function— u(z). The Fourier transform of the bulk 
magnetization function— / (Eq. (45)) 

+ / 3 /jt +I 2 j2k +lj2k 

f=Zi 2L Z %nun^n§\?- n iKy- n 2h^- n &~\ W 

n 3 = 0 »2=-/ 2 /2£ n x =~lj2k 

is 

15 ^=E Z Z X^e-*"^^ (50) 

where x n = njc^ y n = n 2 k 2 ; z n = -ri 3 k 3 . The Fourier transform of 
u(z) = l forz > 0 and u(z) = 1 forz< 0 is [8] 

u{K)=\s{K)^^ (5i) 

The Fourier transform of the system function (See APPENDIX V) 
20 *>> + / + z 2 f~(p 2 +*T (52) 



is 



P * 



aft * t L *^^1 (54) 



25 where 
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Jr 


— 9 irf 


— JL — 






X 


K 


= 2nf y 






y 


K 


= 2nf 2 


= 2k- 

z 



(55) 



The resultant function of the product of the functions H (Eq. (52)) and F (Eq. (49)) then 
convolved with the function U(k^) (Eq. (50)) is 
S = HF® U[k 2 ]= 

(56) 



^Lj ^"1,"2,«3^ 



+ l®iiy y 2 y y ^w.+v.] 



S ^ K z z n =0 y n —l 2 \2 x n =-lj2 

5 The result is given as follows, and the derivation appears in APPENDIX VI. 
S= l -HF+± + f + f X m e< k ' x ^ k 4^ 1 



+ / 2 /2 +/j 



z„ = 0 .y„ =- / 2 / 2 x„ =- l\l 2k 



i' 2 r if 77-fr 2 , i 



7^ 



1 



(57) 



Substitution of the Fourier transform of the system function, H (Eq. (52)), the Fourier 
transform of the bulk magnetization function, F (Eq. (49)), and factoring out /c p 2 in the 
1 0 second term gives 



f2^ 2 2^1^ ^ + " 2 -,[M^V„] 

_z n =Q y n = -h/2 x n =~l\l2 



r 



-/ 3 +/ 2 /2 +/!/2 



2«=0 ^=-/ 2 /2 x n =-hl2k I rJ 1 , -A^ 



r *i 



(58) 

Multiplication by the complex conjugates gives 
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^ 9 ^-J 5j ^'il»«2>K3^ 



-/ 3 +/ a /2 +A/2 



1 


1-7— 


F i -A 1 


1 -A 



r *i 

i+Jr 



+■ 



r *i 



L 1_ ^JL 1+/ ?J 



(59) 



Factoring out the common terms gives 



r 



1+7'- 2 - 



+ 



1 

I" J" 



1 



(60) 



5 Combining exponential terms gives 



1 



(61) 



The function 5 divided by the function H is 



£ j -fe +/ 2 /2 +/,/2 

Tt "™ <7 2 S 2 ^/ii,n 2 ,n 3 



X 



z. <0 



(62) 



10 The inverse Fourier transform of the function S divided by the function # is given as 
follows, where the symbol F _1 (Q) is defined as the inverse Fourier transform of the 
function O. 

-j\k x x„ + k yil +k z 2 n ] I 

I 

}; z„<0 (63) 
The inverse Fourier transform of the function F is 

15 F X (F)=/^X X J! Xn^nAx-^y-yn'Z-Zn] (<*) 



( j ~/ 3 +/ 2 /2 +/,/2 

F -ir_^l =F -^ 8z »=° '»=-'*' 2 *» = - ,l,2 
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Specifically, the inverse Fourier transform of the function F in terms of the dipole 
spacing is 

(~h +y2 +4/2 1 

HE E E % ni ,^3 e " 4 ' ] 

U=0 y n = -hl2 Xn = ~h/2 J 

-/ 3 +/ 2 /2A: +/i/2Jt 

= 2 S X ^ I ,« 2 ,«3 5 [ x ~ n i^-> ; " 7z 2^ z " w 3^] 

ZfJ =0 y„=-t 2 /2k x n = -l\llk 

Consider the general Fourier transform pair in cylindrical coordinates: 
5 g(r,e,z)<^G(s,<j>,(D)(66) 

Under circular symmetry, that is when g is independent of 0 (and hence G independent 
of <f) ) the inverse Fourier transform is [9] 

00 OO 

h{r,z) = 2it\ \H(s,coy 0 (^Ksr)e i2ncoz sdsdco (67) 

0 -co 

The inverse Fourier transform of Eq. (62): 

1 i V V V v -4***" + V« +M- J 1 



Z Z Z ^"l >*2 >«3 



10 " !>; z„<0 (68) 



given by Eq. (64) and Eq. (66) is 



j 



-co 0 

OO 00 



+2* J J e +k ^J 0 \ p p\k p k p ^dk z 



-oo 0 
oo co 



-oo 0 

GO CO 



\-J2tv j je- k ^J 0 fap} p dk 2 e^dk 2 \ 



15 



(69) 

The solution of Eq. (68) appears as follows, and the derivation appears in APPENDIX 
VII. 
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p-l jl ^ ^ \£ x e -j[ k * X n + k yy tt +k **n ] e - k pl*«\ + e +k pl*«\ + j A. H^P^I _ ( 

[8 Z;j = 0 y„=-hl2 x n =-hl2 1,2,2 L 



Z « 3 = 0 n 2 =-/ 2 /2£ m=-hi2k |_ 



Inverse Transform 1 



(70) 



5 Thus, the inverse Fourier transform of — (Eq. (62)) is given by Eq. (69). 

H 



F IhI 



F |oE Z Z ^«i>«2^ 

L°^=o y„=-h/2 x„ = -li/2 



-j[k x x„+k y y n ±k z z n 



+/ 2 /2Ar +/ 1 /2A 



= fZ Z Z %n^n^V-^y" n 2 k ^-^\® 



n 3 = 0 n 2 ~~hi2k «i=-/j/2A 



10 The convolution replaces each coordinate with a spatially shifted coordinate. 

i[sl n r z 1 

2 * *. it- + [ jc - jf -r + D / -^Tj J 

The result of the evaluation of Eq. (71) at each coordinate x n ,y n ,z n is 



15 



(71) 



(72) 



F IhI 



2zl\ 



%n i,«2)«3 

Solving for the bulk magnetization of each dipole gives 



(73) 



F Vh\ 



x „>y„ 



% 
2z 2 



(74) 



The solution of the bulk magnetization of each dipole follows from Eq. (73). 
Discrete values of the voltages produced at the detector array due to the secondary 
magnetic field are recorded during a scan which represent discrete values of function s 
(Eqs. (41-42)); thus, in the reconstruction algorithm that follows, discrete Fourier and 
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Inverse Fourier transforms replace the corresponding continuous functions of Eq. (73) of 
the previous analysis. 

Discrete values of H of Eqs. (52-54), the Fourier transform of the system 
function, replace the values of the continuous function. Furthermore, each sample 
5 voltage of an actual scan is not truly a point sample, but is equivalent to a sample and 
hold which is obtained by inverting the grid antennae matrices or which is read directly 
from a micro-antennae as described in the Finite Detector Length Section. The spectrum 
of a function discretely recorded as values, each of which is equivalent to a sample and 
hold, can be converted to the spectrum of the function discretely recorded as point 
10 samples by dividing the former spectrum by an appropriate sine function. This operation 
is performed and is described in detail in the Finite Detector Length Section. From these 
calculated point samples, the bulk magnetization function (NMR image) is obtained 
following the operations of Eq. (73) as given below. 

15 Reconstruction Algorithm 

Record the RF NMR signal at discrete points over the sample space. Each point is 
designated (x, y, z, RF) and each RF value is an element in matrix A. 

Invert the detector grid matrices defined by the noncommon areas of the overlapping 
elements of the detector array described in the Finite Detector Length Section to obtain 
the sample and hold voltages which form Matrix A 1 (if micro-antennas are used, form 
Matrix A ? of the recorded voltages directly). 

Fourier transform the time dependent signals to give the intensity and phase of each 
component. The NMR signal of each voxel at any given detector gives rise to a Fourier 
component with a unique phase angle relative to the Fourier component of any other 
voxel of the phantom at that detector. The matrix of Fourier components that correspond 
to the NMR signal of a given voxel over the detectors is determined. This may be 
achieved by using a first component having a phase angle and calculating the phase angle 
as a function of spatial position of the first detector relative to any other detector and 
identifying the component at each detector having the calculated phase angle. The 
matrices A n are determined for all of the voxels. The measurements of the spatial 
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variations of the transverse RF field of a given matrix is used to determine the coordinate 
location of each voxel. Thus, each matrix of components associated by phase comprises 
the intensity variation over the sample space of the RF field of the bulk magnetization M 
of each voxel. 

5 

4) Three-dimensionally Fourier transform each matrix A n , using a discrete Fourier 
transform formula such as that which appears in McC. Siebert [10] to obtain each matrix 
B n of elements at frequencies corresponding to the spatial recording interval in each 
direction. 

10 A n (x,y,z)=>B n (k x ,k y ,k z ) 

5) Multiply each element of each matrix B n by a value which is the inverse of the Fourier 
transform of a square wave evaluated at the same frequency as the element where the 
square wave corresponds to a sample and hold operation performed on the continuous 

15 voltage function produced at the detector array by the secondary field. This 

multiplication forms each matrix B n * which is the discrete spectrum of the continuous 
voltage function discretely point sampled (See the Finite Detector Length Section for 
details of this operation). 

20 

6) Multiply each element of each matrix B n * by the value which is the inverse (reciprocal) 
of the Fourier transform of the system function evaluated at the same frequency as the 
element to form each matrix B n * * . 

25 

7) Inverse three-dimensionally Fourier transform each matrix B n * * using the discrete 
inverse Fourier transform formula such as that which appears in McC. Siebert [10] to 
form each matrix C n whose elements correspond to the bulk magnetization of the dipoles 
at the points of the image space of spatial interval appropriate for the frequency spacing 

30 of points of matrix B n * * . 

B n **(k x ,k y ,k 2 )=>C n (x,y,z) 
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7F 

8) Divide each element of each matrix C n by — - to correct for the restriction that the 

2Z n 

sample space is defined as z greater than zero (z> 0). This operation creates each matrix 
D n which is the bulk magnetization M map (NMR image). 



2z 



2 



In the NMR case shown in FIGURES 8, 10, 11, and 12, Step 7 of the 
Reconstruction Algorithm of the Fourier Transform Reconstruction Algorithm Section is 
10 equivalent to the general case except for the coordinate designations and the relationship 
between the magnetic moment of a static secondary field and the bulk magnetization M 
of the RF field wherein m z , the magnetic moment along the z-axis, of Eq. (1) 
corresponds to the bulk magnetization M of each voxel: 



8*1)5 Divide each element of each matrix C n by — - to correct for the restriction that the 

sample space is defined as y greater than zero, y > 0 . This operation creates each matrix 
D n which is the bulk magnetization M map. 



In one embodiment, the point spread of the reconstructed voxel is corrected by assigning 
20 one voxel above a certain threshold with the bulk magnetization M . The other voxels 
are assigned a zero value. This procedure may be repeated for all voxels. In the limit 
with sufficient phase resolution, each volume element is reconstructed independently in 
parallel with all other volume elements such that the scan time is no greater than the 
nuclear free induction decay (FID) time. 

25 

9) Superimpose the separate bulk magnetization M maps corresponding to each matrix A n 
Plot the bulk magnetization M with the same spatial interval in each direction as the 
sampling interval in each direction over all matrices A n (i.e.. the total bulk magnetization 
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M map which is plotted and displayed is the superposition of the separate maps of each 
magnetic moment corresponding a matrix A n ). 

(The above steps for each matrix relate generally to the program implementation shown 
5 in the listing of the Exemplary Reconstruction Program as follows. The above Steps 1 
and 2 relate to the Data Statements beginning at lines 50; and Step 4 relates to the X, Z 
and Y FFT operations of lines 254, 455 and 972, respectively. Steps 5 and 6 are 
implemented by the processes of lines 2050, 2155 and 2340; and Step 7 relates to the X, 
Y and Z inverse transform of lines 3170, 3422 and 3580, respectively. Step 8 and 8* . 
10 relates to the correction and normalization process of line 4677.) 



Exemplary Reconstruction Program 

10 ! 4D-MRI ALGORITHM 

25 OPEN "4D-MRI PROTOTYPE I . LIS" FOR OUTPUT AS #1 
15 26 C= . 05 

27 DIM X(9,9,9) 

2 8 FOR W=l TO 9 STEP 1 

2 9 FOR T=l TO 9 STEP 1 

3 0 FOR S = l TO 9 STEP 1 
20 31 X(W,T,S)=0 

32 NEXT S 
3 3 NEXT T 

34 NEXT W 

35 X(5,5, 6) =1 
25 36 DIM DI (9,9) 

3 7 PRINT #1, "DIPOLE PHANTOM" 

40 FOR Q=l TO 9 STEP 1 

41 FOR R=l TO 9 STEP 1 

42 FOR U=l TO 9 STEP 1 
30 43 LET H=X(U,R,Q) 

44 LET DI(U,R)=H 

45 NEXT U 

4 6 NEXT R 

47 MAT PRINT #1,DI, 
35 4 8 NEXT Q 

5 0 DATA 0,0,0,0,0,0,0,0,0,0,0,0,0 

52 DATA 0,0,0,0,0, -.6, -.9, -1.1, -1.2,-1.8,5.7,3.5, .4, -.4, -.5, -.7, -.6, .7 
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54 DATA 4. 2, 19. 7,7.8, .1,-. 9, -.3, -.4, -.7, .7, 13. 9, 25. 65, 10.1, .9,-. 8,-. 3,- 
.6 

56 DATA -.4, 2. 7, 10. 7, 12. 4, 5. 6, 2. 6, -.5, -.3, -.5, -.6, -.6, .3, 3. 7, .9, -.5,- 
1.0 
5 58 DATA 

0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0 
60 DATA 0 , 0, - . 5, - . 8 , - . 8, -1.0 , -1.2, - .8, 6 .9, .2, - .2, - .4, - . 6, - 
.6,1.2,4.1,12.7,0, .6 

62 DATA -. 6, -.3,-. 5, -.7, -.3, 3. 7, 9. 05, 7. 6, 1.3, -.5,-. 3,-. 4,- 
10 .2,2.7,4.3,10.6, .8 

64 DATA 1.7, -.4, -.3, -.5, -.4, -.1,1. 1,1. 1,1. 7,. 2,- 
.6,0,0,0,0,0,0,0,0,0,0,0,0,0 

66 DATA 0, 0,0, 0,0, 0,0, 0,0, 0,0, 0,0, 0,0, 0,0, 0,0,0,0, 0,0, -.5,-. 6,-. 6,-. 5 , - 
.5 

15 68 DATA .3, 2. 5, .l,-.l,-. 3, -.4, -.3, 1.3, 3. 9, 8. 1,1. 7, .9, -.3, -.3, -.4,-. 6,- 
.8,1.5 

70 DATA 1.65, 4. 7, 1.3, -.3, -.2, -.2, .1,2. 1,3. 5, 7.1, .7, .7, -.3, -.2, -.4,- 
.2, .2,1.2 

72 DATA 1.1,1.4, . 4, - 
20 .3,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0 

74 DATA 0, 0,0, 0,0, 0,0, 0,0,0, -.4, -.5,0,-.!, -.1,1. 7, 1.5,. 2,-. 1,-. 2, - 
.2,0,1.3,3.1 

76 DATA 4. 1, 3. 4,. 9, 0, 2, -.3, -.4, -.4, 1.1,. 65, 2. 9, 1.1,-. 2, -.1,- 
.1, .2,1.2,1.6 

25 78 DATA 4.9, .1, .5, -.3, -.2,-. 3, -.1,0, .3,1.3, .6, .1, - 
.1,0,0,0,0,0,0,0,0,0,0,0,0 

8 0 DATA 0, 0,0, 0,0, 0,0, 0,0, 0,0, 0,0, 0,0, 0,0, 0,0,0, 0,0,0, 0,-. 3 , - 
.2, .1, .9,1.6,1.6 

82 DATA . 4,0,0, -.l,-.l,. 2, 1.0, 2. 2, 2. 7, 3. 7,. 5,0, -.2,-. 2,-. 2,- 
30 .1, .6, .25,1.6, .7 

84 DATA - .2, .1, .7, .9,2.8,0, .3,-.2,-.l,-.2,-.l,0, .2, .3, .3,0,- 
.1,0,0,0 

85 DATA 0,0,0,0,0,0,0,0,0,0,0,0,0,0,0 

86 DATA 0, 0,0, 0,0, 0,0, 0,0, 0,0, 0,0, 0,0, 0,0,0, -.2,- 
35 .1,. 1,1. 0,1. 4, 1.3,. 2, 0,0 

88 DATA - . 1, 0, .2 , . 7, 1.5,2 . 0, 2 . 9, . 3, 0, - . 1, - . 1, . 1, 0, . 7, .25, .9, .4 , - 
.1,0,0, .1, .5 

90 DATA .5, 1-3,0, .2, -.1,0, -.1,0, .1,-2, .1,-2,0, - 
.1,0,0,0,0,0,0,0,0,0,0,0,0,0,0 
40 92 DATA 0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,- 
.1,0, .1, .8, .9,1.0, .1,0,0 
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94 DATA 0,0, .2, .5, 1.1, 1.3,2.0, .3, 0,0, -.1,0,0, .3, .15, .6, .2,- 
.1,0,0, .1, .4, .3, .9 
96 DATA 0, . 2, -.1, 0, - 

.1, .1, .1, .1, .1, .2,0,0,0; 0,0,0,0,0,0,0,0,0,0,0,0,0,0,0 
5 98 DATA 0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,- 
.1,0,0, .7, .7, .6, .1,0, .1 

100 DATA 0, .1,. 2, .3, .6, .8, 1.0, .2,0,0,- 
.1,0,0, .1, .15, .3, .1,0,0,0,0, .2, .2, .5 
102 DATA 0, . 1, - 

10 .1,0,0,0, .1, .1, .1, .1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0 
104 DATA 

0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0, .5, .4, .4, .1,0,0,0,0, .1, .1 
106 DATA .4, .5, .6, .1,0,0,0,0,0,0, .05, .1,0,0,0,0,0, .1, .1, .3,0, . 1, - 
.1,0,0,0,0,0 

15 108 DATA .1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0 
120 DIM V(9, 9, 9) 

13 0 FOR Z=l TO 9 STEP 1 

14 0 FOR Q=l TO 9 STEP 1 
145 FOR U=l TO 9 STEP 1 

20 147 READ V(Z,Q,U) 

222 NEXT U 

223 NEXT Q 

224 NEXT Z 

22 5 PRINT #1, "VOLTAGE DATA" 
25 226 DIM VO (9, 9) 

227 FOR Q=l TO 9 STEP 1 
22 8 FOR R=l TO 9 STEP 1 

229 FOR U=l TO 9 STEP 1 

230 LET H=V(U,R,Q) 
30 231 LET VO(U,R)=H 

232 NEXT U 

233 NEXT R 

234 MAT PRINT #1, VQ , 

235 NEXT Q 

35 254 ! FFT THE ROWS OF SAMPLED VOLTAGES IN THE X DIRECTION 

255 DIM MR(9) 

256 DIM MI (9) 

257 DIM R(9) 

258 DIM RV(9,9,9) 
40 259 DIM IV(9,9,9) 

260 FOR V=0 TO 9 STEP 1 
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270 FOR M=l TO 9 STEP 1 
2 80 FOR N=l TO 9 STEP 1 

2 90 LET Y=V(M,N,V) 
300 LET R (N) =Y 

5 3 05 NEXT N 

310 EXTERNAL SUB FFT (DIM ( ) , DIM ( ) , DIM ( ) ) 
32 0 CALL FFT (R { ) , MR ( ) , MI ( ) ) 
330 FOR N=l TO 9 STEP 1 
34 0 LET Y=MR(N) 
10 350 LET RV(M,N,V)=Y 

3 60 NEXT N 

370 FOR N=l TO 9 STEP 1 

3 80 LET G=MI (N) 
'390 LET IV(M,N,V)=G 
15 4 00 NEXT N 

410 NEXT M 

415 NEXT V 

417 DIM RVO (9, 9) 

41*9 DIM IVO(9,9) 
20 420 PRINT #1,"RV" 

421 FOR Q=l TO 9 STEP 1 

422 FOR R=l TO 9 STEP 1 

423 FOR U=l TO 9 STEP 1 

424 LET H=RV(U / R / Q) 
25. 425 LET RVO(U,R)=H 

42 6 NEXT U 
42 7 NEXT R 

42 8 MAT PRINT #1, RVO, 

42 9 NEXT Q 

30 43 0 PRINT #1, "IV" 

431 FOR Q=l TO 9 STEP 1 

432 FOR R=l TO 9 STEP 1 

433 FOR U=l TO 9 STEP 1 

434 LET H=IV(U,R,Q) 
35 43 5 LET IVO(U,R)=H 

436 NEXT XJ 

43 7 NEXT R 

43 8 MAT PRINT #1, IVO, 
439 NEXT Q 

40 45 5 I FFT THE COLUMNS OF THE SAMPLED VOLTAGE IN THE Z DIRECTION 
460 DIM RRV(9, 9, 9) 
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470 


DIM IRV(9,9,9) 






475 


FOR V= 1 TO 9 STEP 


1 




480 


FOR N=l TO 9 STEP 


1 




490 


FOR M=l TO 9 STEP 


1 


5 


500 


LET Y=RV(M,N,V) 






510 


LET R (M) =Y 






520 


NEXT M 






530 


CALL FFT(R() ,MR() ,MI() ) 




540 


FOR M=l TO 9 STEP 


1 


10 


550 


LET H=MR (M) 






560 


LET RRV(M,N,V)=H 






570 


NEXT M 






580 


FOR M=l TO 9 STEP 


1 




590 


LET G=MI (M) 




15 


600 


LET IRV(M,N,V)=G 






610 


NEXT M 






620 


NEXT N 






625 


NEXT V 






626 


DIM RRVO (9,9) 




20 


627 


DIM IRVO(9,9) 






630 


PRINT #1, " RRV " 






631 


FOR Q=l TO 9 STEP 


1 




632 


FOR R=l TO 9 STEP 


1 




633 


FOR U =1 TO 9 STEP 1 


25 


634 


LET H=RRV(U,R,Q) 






637 


LET RRVO (U,R) =H 






638 


NEXT U 






639 


NEXT R 






640 


MAT PRINT #1, RRVO, 


30 


641 


NEXT Q 






650 


PRINT #1, "IRV" 






651 


FOR Q=l TO 9 STEP 


1 




653 


FOR R=l TO 9 STEP 


1 




654 


FOR U=l TO 9 STEP 


1 


35 


655 


LET H= IRV(U,R,Q) 






656 


LET IRVO(U,R)=H 






657 


NEXT U 






658 


NEXT R 






659 


MAT PRINT #1, IRVO, 


40 


660 


NEXT Q 






661 


DIM RIV(9,9,9) 
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662 DIM IIV(9,9,9) 

685 FOR V=l TO 9 STEP 1 

690 FOR N=l TO 9 STEP 1 

700 FOR M=l TO 9 STEP 1 

5 710 LET Y=IV(M,N,V) 

760 LET R(M)=Y 

770 NEXT M 

830 CALL FFT (R ( ) , MR ( ) , MI ( ) ) 

84 0 FOR M=l TO 9 STEP 1 

10 85 0 LET H=MR(M) 

860 LET RIV(M,N, V) =H 

87 0 NEXT M 

8 72 FOR M= 1 TO 9 STEP 1 

87 3 LET G=MI (M) 

15 874 LET IIV(M,N,V)=G 

875 NEXT M 

876 NEXT N 
8 77 NEXT V 

878 DIM RIVO(9, 9) 

20 879 DIM IIVO(9,9) 

88 0 PRINT #1, "RIV" 

8 81 FOR Q=l TO 9 STEP 1 

8 82 FOR R=l TO 9 STEP 1 

8 83 FOR U=l TO 9 STEP 1 

25 884 LET H=RIV(U,R,Q) 

8 85 LET RIVO(U,R)=H 

886 NEXT U 

887 NEXT R 

8 88 MAT PRINT #1, RIVO, 

30 889 NEXT Q 

890 PRINT #1, "IIV" 

891 FOR Q=l TO 9 STEP 1 
8 92 FOR R=l TO 9 STEP 1 
893 FOR U=l TO 9 STEP 1 

35 894 LET H=IIV(U,R,Q) 

895 LET IIVO(U,R)=H 

896 NEXT U 
8 97 NEXT R 

898 MAT PRINT #1, IIVO, 

40 899 NEXT Q 

900 DIM RVA(9, 9,9) 
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901 


DIM IVA(9,9,9) 






904 


DIM RVAO(9,9) 






906 


DIM IVAO(9,9) 






908 


FOR Q=l TO 9 STEP 


1 


5 


910 


FOR R=l TO 9 STEP 


1 




911 


FOR U= 1 TO 9 STEP 


1 




912 


LET H=IIV(U,R,Q) 






913 


H= (-1) *H 






914 


LET G=RRV(U / R / Q) 




10 


915 


LET L-G+H 






916 


LET RVA(U,R,Q) =L 






917 


NEXT U 






918 


NEXT R 






919 


NEXT Q 




15 


920 


FOR Q-l TO 9 STEP 


1 




921 


FOR R=l TO 9 STEP 


1 




922 


FOR U=l TO 9 STEP 


1 




923 


LET H=IRV(U,R,Q) 






924 


LET L=RIV(U,R,Q) 




20 


925 


LET IVA(U,R,Q) =H+L 




927 


NEXT U 






928 


NEXT R 






930 


NEXT Q 






932 


PRINT #1, "RVA" 




25 


934 


FOR Q=l TO 9 STEP 


1 




936 


FOR R=l TO 9 STEP 


1 




938 


FOR U=l TO 9 STEP 


1 




940 


LET H=RVA(U,R,Q) 






942 


LET RVAO(U,R)=H 




30 


944 


NEXT U 






946 


NEXT R 






948 


MAT PRINT #1, RVAO, 




950 


NEXT Q 






952 


PRINT #1, "IVA" 




35 


954 


FOR Q=l TO 9 STEP 


1 




956 


FOR R=l TO 9 STEP 


1 




958 


FOR U= 1 TO 9 STEP 


1 




960 


LET H=IVA(U / R / Q) 






962 


LET IVAO(U,R)=H 




40 


964 


NEXT U 






966 


NEXT R 
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66 

968 MAT PRINT #1, IVAO, 

97 0 NEXT Q 

972 1FFT THE ROWS OF THE SAMPLED VOLTAGES IN THE Y DIRECTION 

98 0 DIM RVAY (9, 9, 9) 
5 990 DIM IRVAY(9,9,9) 

992 DIM RVAYO(9,9) 
994 DIM IRVAYO(9,9) 

1012 FOR M=l TO 9 STEP 1 

1013 FOR N=l TO 9 STEP 1 
10 1014 FOR V=l TO 9 STEP 1 

1015 LET Y=RVA(M,N,V) 

1016 LET R(V) =Y 

1017 NEXT V 

1018 CALL FFT(R() , MR ( ) , MI () ) 
15 1019 FOR V=l TO 9 STEP 1 

102 0 LET H=MR(V) 

1021 LET RVAY (M,N,V)=H 

1022 NEXT V 

1023 FOR V=l TO 9 STEP 1 
20 1024 LET G=MI(V) 

1025 LET IRVAY(M / N / V) =G 

102 6 NEXT V 

102 7 NEXT N 

102 8 NEXT M 
25 1030 PRINT #1,"RVAY" 

1032 FOR Q=l TO 9 STEP 1 

1034 FOR R=l TO 9 STEP 1 

1036 FOR U=l TO 9 STEP 1 

1038 LET H=RVAY(U,R,Q) 
30 1040 LET RVAYO (U,R) =H 

1042 NEXT U 

1044 NEXT R 

1046 MAT PRINT #l,RVAYO, 

104 8 NEXT Q 
35 105 0 PRINT #1, " I RVAY" 

1052 FOR Q=l TO 9 STEP 1 

1054 FOR R=l TO 9 STEP 1 

1056 FOR U=l TO 9 STEP 1 

1058 LET H=IRVAY(U,R,Q) 
40 1060 LET IRVAYO (U,R) =H 

1062 NEXT U 
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1064 NEXT R 

1066 MAT PRINT #l,IRVAYO, 
1068 NEXT Q 
1070 DIM RIVY(9,9,9) 
5 1080 DIM IIVY(9,9,9) 

1085 DIM RIVYO(9,9) 

1086 DIM IIVYO(9,9) 
1090 FOR M=l TO 9 STEP 1 
1100 FOR N=l TO 9 STEP 1 

10 113 8 FOR V=l TO 9 STEP 1 
1139 LET Y=IVA(M / N / V) 

114 0 LET R(V) =Y 

1141 NEXT V 

1142 CALL FFT (R ( ) , MR ( ) f MI ( ) ) 
15 1143 FOR V=l TO 9 STEP 1 

1144 LET H=MR(V) 

1145 LET RIVY (M, N, V) =H 

1146 NEXT V 

1147 FOR V=l TO 9 STEP 1 
20 114 8 LET G=MI (V) 

1149 LET IIVY(M,N,V)=G 

115 0 NEXT V 

1151 NEXT N 

1152 NEXT M 

25 1153 PRINT #1, "RIVY" 

1160 FOR Q=l TO 9 STEP 1 

1162 FOR R=l TO 9 STEP 1 

1164 FOR U=l TO 9 STEP 1 

1166 LET H=RIVY(U / R / Q) 
30 1170 LET RIVYO(U,R)=H 

1172 NEXT U 

1174 NEXT R 

117 8 MAT PRINT #1, RIVYO, 

1180 NEXT Q 
35 1185 PRINT #1, "IIVY" 

1190 FOR Q=l TO 9 STEP 1 

12 0 0 FOR R=l TO 9 STEP 1 

1210 FOR U=l TO 9 STEP 1 

1212 LET H= IIVY(U / R / Q) 
40 1214 LET IIVYO(U,R)=H 

1216 NEXT U 
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1218 NEXT R 

122 0 MAT PRINT #l,IIVYO, 

1222 NEXT Q 

1230 DIM YRVA(9 ,9,9) 
5 1240 DIM YIVA(9,9,9) 

1241 FOR Q=l TO 9 STEP 1 

1250 FOR R=l TO 9 STEP 1 

1260 FOR U=l TO 9 STEP 1 

1270 LET L=IIVY(U,R,Q) 
10 1280 LET B=RVAY(U,R,Q) 

1290 LET YRVA(U,R,Q) =B-L 

13 00 NEXT U 

1310 NEXT R 

132 0 NEXT Q 
15 1330 FOR Q=l TO 9 STEP 1 

1340 FOR R=l TO 9 STEP 1 

1345 FOR U=l TO 9 STEP 1 

1350 LET H=RIVY(U,R,Q) 

1360 LET L=IRVAY(U,R,Q) 
20 13 7 0 LET YIVA(U,R,Q)=L+H 

13 8 0 NEXT U 

13 90 NEXT R 

14 00 NEXT Q 

1410 PRINT #1, " YRVA" 
25 1412 DIM YRVAO(9,9) 
1414 DIM YIVAO(9,9) 

142 0 FOR Q=l TO 9 STEP 1 

143 0 FOR R=l TO 9 STEP 1 
1440 FOR U=l TO 9 STEP 1 

30 1450 LET H=YRVA(U,R,Q) 
14 60 LET YRVAO (U,R) =H 
14 70 NEXT U 
14 80 NEXT R 

1490 MAT PRINT #1, YRVAO, 
35 15 0 0 NEXT Q 

1510 PRINT #1, "YIVA" 

152 0 FOR Q=l TO 9 STEP 1 

153 0 FOR R=l TO 9 STEP 1 
1540 FOR U=l TO 9 STEP 1 

40 1545 LET H=YIVA(U,R,Q) 
1550 LET YIVAO(U,R)=H 
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1560 NEXT U 
157 0 NEXT R 

15 80 MAT PRINT #l,YIVAO, 
1590 NEXT Q 

5 2050 1 GENERATE THE DISCRETE SPECTRUM OF THE SYSTEM FUNCTION AND THE 
SINC 

2055 L FUNCTION OF THE SAMPLE AND HOLD CORRESPONDING TO THE FINITE 
DETECTOR 
2 057 ! DIMENSIONS 
10 2060 DIM SFH(9,9,9) 

2061 DIM SFHO(9 / 9) 

2062 DIM SINC (9,9,9) 

2063 DIM SINCO(9,9) 
2065 LET SFH(4 f 4,4) =4*PI 

15 2 066 FOR M=-4 TO 4 STEP 1 
2070 FOR V=-4 TO 4 STEP 1 

2080 FOR N=-4 TO 4 STEP 1 

2081 H=.04 

2082 J=ABS (M) +ABS (N) +ABS (V) 
20 2 085 IF J=0 THEN GO TO 2098 

2 090 T=4*PI* ( (2*PI*N/9*1/ . 1) A 2+ ( (2*PI*V/9*l/ . 1) A 2) ) 
2 092 B= (2*PI*N/9*l/ . 1) ^2+ (2*PI*M/9*l/C) ^2+ (2*PI*V/9*l/ . 1) A 2 
2095 LET SFH(M+5,N+5,V+5) =T/B 
2 098 G=ABS (N) +ABS (V) 
25 2100 IF G=0 THEN GO TO 2114 
2101 A=l 

2105 IF N=0 THEN GO TO 2107 

2106 A=SIN(2*PI*10*N/9*.l)/ (PI*N/9*10) 

2107 B=l 

30 2108 IF V=0 THEN GO TO 2110 

2109 B=SIN(2*PI*10*V/9*.l)/ (PI*V/9*10) 

2110 H=A*B 

2111 IF N=0 THEN H=.2*H 

2112 IF V=0 THEN H=.2*H 

35 2114 LET SINC (M+5,N+5, V+5) =H 
213 0 NEXT N 

2131 NEXT V 

2132 NEXT M 

2135 PRINT #1, "SFH " 
40 2136 FOR Q=l TO 9 STEP 1 
2137 FOR R=l TO 9 STEP 1 
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213 8 FOR U=l TO 9 STEP 1 

2139 LET S=SFH(U,R,Q) 

2140 LET SFHO(U,R)=S 

2141 NEXT U 
5 2142 NEXT R 

2143 MAT PRINT #1, SFHO, 

2144 NEXT Q 

2145 PRINT #1, "SINC" 

214 6 FOR Q=l TO 9 STEP 1 
10 2147 FOR R=l TO 9 STEP 1 

2148 FOR U=l TO 9 STEP 1 

2149 LET S=SINC(U,R,Q) 

2150 LET SINCO(U,R)=S 

2151 NEXT U 
15 2152 NEXT R 

2153 MAT PRINT #1, SINCO, 

2154 NEXT Q 

2155 I INVERSE THE DISCRETE SPECTRUM OF THE SYSTEM FUNCTION AND THE SINC 
2157 ! FUNCTION 

20 2160 DIM HR (9,9,9) 

2170 DIM HRO (9, 9) 

2171 DIM SINCR(9,9,9) 
2175 FOR V=l TO 9 STEP 1 
2180 FOR M=l TO 9 STEP 1 

25 2190 FOR N=l TO 9 STEP 1 

2200 LET Y=SFH(M,N,V) 

2201 LET H=SINC (M, N, V) 

2210 IF Y=0 THEN GO TO 2212 

2211 Y=l/Y 

30 2212 IF H=0 THEN GO TO 2230 

2221 H=l/H 

2230 LET HR(M,N,V)=Y 

2235 LET SINCR (M, N, V) =H 

2 24 0 NEXT N 

35 2250 NEXT M 

2260 NEXT V 

2310 PRINT #1, "HR" 

2 311 FOR Q=l TO 9 STEP 1 

2 312 FOR R=l TO 9 STEP 1 

40 2 313 FOR U=l TO 9 STEP 1 

2314 LET S=HR(U,R,Q) 
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2 315 LET HRO(U,R)=S 
2 316 NEXT U 
2317 NEXT R 
2 32 0 MAT PRINT #l,HRO, 
5 2321 NEXT Q 

2340 ! DIVIDE THE TRANSFORMED DATA BY THE TRANSFORM OF THE SYSTEM 





FUNCTION 






2345 


1 AND THE SINC FUNCTION 




3030 


DIM FR (9, 9, 9) 




10 


3050 


DIM FI (9, 9, 9) 






3052 


DIM FRO (9,9) 






3054 


DIM FIO (9,9) 






3056 


FOR V=l TO 9 STEP 


1 




3057 


FOR M=l TO 9 STEP 


1 


15 


3058 


FOR N=l TO 9 STEP 


1 




3059 


T=YRVA(M / N / V) 






3060 


S=HR(M,N, V) 






3062 


L=SINCR(M,N, V) 






3065 


K=S*T*L 




20 


3066 


LET FR(M / N,V)=K 






3067 


NEXT N 






3068 


NEXT M 






3070 


NEXT V 






3080 


FOR V=l TO 9 STEP 


1 


25 


3081 


FOR M=l TO 9 STEP 


1 




3082 


FOR N=l TO 9 STEP 


1 




3083 


H=SINCR (M,N, V) 






3093 


K= HR(M,N,V) 






3094 


L= YIVA(M / N / V) 




30 


3095 


E=K*L*H 






3096 


LET FI (M,N, V) =E 






3097 


NEXT N 






3098 


NEXT M 






3100 


NEXT V 




35 


3130 


PRINT #1, "FR" 






3131 


FOR Q=l TO 9 STEP 


1 




3132 


FOR R=l TO 9 STEP 


1 




3133 


FOR U=l TO 9 STEP 


1 




3134 


LET H= FR(U,R,Q) 




40 


3135 


LET FRO(U,R)=H 






3136 


NEXT U 
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313 7 NEXT R 

313 8 MAT PRINT #1, FRO, 

3140 NEXT Q 

3141 PRINT #1, " FI 11 

5 3142 FOR Q=l TO 9 STEP 1 

3143 FOR R=l TO 9 STEP 1 

3144 FOR U=l TO 9 STEP 1 

3145 LET H=FI(U,R,Q) 

3146 LET FIO(U,R)=H 
10 3156 NEXT U 

3157 NEXT R 

315 8 MAT PRINT #1, FIO, 
3160 NEXT Q 

317 0 1 INVERSE TRANSFORM THE ROWS IN THE X DIRECTION 
15 3180 DIM RF(9,9,9) 

3185 DIM RFO(9,9) 

3187 DIM MIFO (9,9) 

3190. DIM MIF(9,9,9) 

3195 DIM 1(9) 
20 3196 FOR V=l TO 9 STEP 1 

3200 FOR M=l TO 9 STEP 1 

3210 FOR N=l TO 9 STEP 1 

3220 LET Y=FR(M,N,V) 

3230 LET R(N) =Y 
25 3240 NEXT N 

3250 FOR N=l TO 9 STEP 1 

3260 LET Y=FI(M,N,V) 

32 7 0 LET I(N)=Y 

3280 NEXT N 

30 32 85 EXTERNAL SUB IFT (DIM ( ) , DIM ( ) , DIM ( ) , DIM ( ) ) 

3290 CALL IFT (R ( ) , I ( ) , MR ( ) , MI ( ) ) 

3300 FOR N=l TO 9 STEP 1 

3 310 LET Y=MR(N) 

3320 LET RF(M,N,V)=Y 
35 3330 NEXT N 

3340 FOR N=l TO 9 STEP 1 

335 0 LET Y=MI (N) 

3360 LET MIF (M,N, V) =Y 

3370 NEXT N 
40 3380 NEXT M 

3385 NEXT V 
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3390 PRINT #1, "RF" 

3391 FOR Q=l TO 9 STEP 1 

3392 FOR R=l TO 9 STEP 1 

3393 FOR U=l TO 9 STEP 1 
5 3394 LET H=RF(U,R,Q) 

3 3 95 LET RFO(U,R)=H 

3396 NEXT U 

33 97 NEXT R 

3400 MAT PRINT #1, RFO, 

10 3405 NEXT Q 

3410 PRINT #1," MIF " 

3411 FOR Q=l TO 9 STEP 1 

3412 FOR R=l TO 9 STEP 1 

3413 FOR U=l TO 9 STEP 1 
15 3414 LET H= MIF (IT, R, Q) 

3415 LET MIFO(U,R)=H 

3416 NEXT U 

3417 NEXT R 

342 0 MAT PRINT #1, MIFO, 

20 3421 NEXT Q 

3422 ! INVERSE TRANSFORM THE ROWS IN THE Y DIRECTION 

3430 DIM RF Y (9,9,9) 

3432 DIM MIFY(9 / 9 / 9) 

3433 DIM RFYO (9, 9) 
25 3434 DIM MIFYO(9,9) 

3435 FOR M=l TO 9 STEP 1 

3440 FOR N==l TO 9 STEP 1 

3450 FOR V=l TO 9 STEP 1 

3460 LET Y=RF(M,N,V) 

30 3470 LET R(V)=Y 

3480 NEXT V 

3490 FOR V=l TO 9 STEP 1 

3500 LET Y=MIF(M,N,V) 

3510 LET I (V) =Y 

35 3 52 0 NEXT V 

3525 EXTERNAL SUB IFT (DIM ( ) , DIM ( ) , DIM ( ) , DIM ( ) ) 

3526 CALL IFT (R ( ) , I ( ) , MR ( ) , MI ( ) ) 

3527 FOR V=l TO 9 STEP 1 
3 52 8 LET Y=MR (V) 

40 3529 LET RFY (M, N, V) =Y 

3530 NEXT V 
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3531 FOR V=l TO 9 STEP 1 

3532 LET Y=MI (V) 

3533 LET MIFY (M,N, V) =Y 

3534 NEXT V 
5 3535 NEXT N 

3536 NEXT M 

3 53 7 PRINT #1, "RFY" 

353 8 FOR Q=l TO 9 STEP 1 
3539 FOR R=l TO 9 STEP 1 

10 3540 FOR U=l TO 9 STEP 1 

3541 LET H=RFY(U,R,Q) 

3542 LET RFYO(U,R)=H 

354 3 NEXT U 
3544 NEXT R 

15 3545 MAT PRINT #1, RFYO , 

3546 NEXT Q 

3547 PRINT #1, "MIFY" 
3550 FOR Q=l TO 9 STEP 1 
3 555 FOR R=l TO 9 STEP 1 

20 3560 FOR U=l TO 9 STEP 1 

3565 LET H=MIFY (U, R, Q) 

3566 LET MIFYO (U, R) =H 

3567 NEXT U 

3568 NEXT R 

25 3570 MAT PRINT #1, MIFYO, 
3 575 NEXT Q 

3 5 8 0! INVERSE TRANSFORM THE COLUMNS IN THE 2 DIRECTION 
3581 DIM F(9,9,9) 
3590 DIM FO(9,9) 
30 3592 FOR V=l TO 9 STEP 1 

3593 FOR N=l TO 9 STEP 1 

3594 FOR M=l TO 9 STEP 1 
3 600 LET Y=RFY (M, N, V) 
4546 LET R (M) =Y 

35 4547 NEXT M 

454 8 FOR M=l TO 9 STEP 1 

4549 LET Y=MIFY (M , N, V) 

4550 LET I (M) =Y 
4552 NEXT M 

40 4553 EXTERNAL SUB IFT2 (DIM ( ) , DIM ( ) , DIM ( ) , DIM ( ) ) 
4554 CALL IFTZ (R ( ) , I ( ) , MR ( ) , MI ( ) ) 
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4556 FOR M=l TO 9 STEP 1 

4557 LET Y=MR (M) 
4559 LET F(M,N,V)=Y 
4600 NEXT M 

5 4602 NEXT N 

4604 NEXT V 

4605 PRINT #1, "F" 

4610 FOR Q=l TO 9 STEP 1 

4 62 0 FOR R=l TO 9 STEP 1 
10 4630 FOR U=l TO 9 STEP 1 

4635 LET H=F(U,R,Q) 

4640 LET FO(U,R)=H 

4 65 0 NEXT U 

4660 NEXT R 
15 4666 MAT PRINT #l,FO, 

4670 NEXT Q 

4677 "CORRECT FOR THE U(Z) CONVOLUTION AND NORMALIZE THE RECONSTRUCTION 

4678 DIM CF (9, 9,9) 
4780 DIM CFO (9, 9) 

20 4809 FOR V=l TO 9 STEP 1 

4 810 FOR N=l TO 9 STEP 1 

4871 LET K=F(1,N,V) 

4873 K=K/10E8 

4974 LET CF(1,N,V)=K 
25 5915 NEXT N 

6616 NEXT V 

6617 FOR V=l TO 9 STEP 1 

662 0 FOR M=2 TO 9 STEP 1 

663 0 FOR N=l TO 9 STEP 1 
30 6640 LET K=F(M,N,V) 

6642 P=0.5*PI/ (M-l) ^2/C^2 
6650 K=K/P 

6660 LET CF(M,N,V)=K. 
667 0 NEXT N 
35 668 0 NEXT M 
6681 NEXT V 

6690 PRINT #1, "RECONSTRUCTION" 

6691 FOR Q=l TO 9 STEP 1 

6692 FOR R=l TO 9 STEP 1 
40 6693 FOR U=l TO 9 STEP 1 

6694 LET H=CF(U,R,Q) 
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6695 H=H/CF(5,5,6) 

6696 LET CFO(U,R)=H 

6697 NEXT U 

6698 NEXT R 

5 67 0 0 MAT PRINT #l,CFO, 
6705 NEXT Q 
6710 END 

6722 SUB PPT(R() , MR ( ) , MI ( ) ) 
672 6 FOR M=l TO 9 STEP 1 
10 6730 A=0 

6740 FOR N=l TO 9 STEP 1 
6750 LET H=R(N) 

6760 B=H*COS (2*PI* (M-5) * (N-5) / 9) 
677 0 A=A+B 
15 6780 NEXT N 
6790 A=A/9 
6800 LET MR (M) =A 
6810 NEXT M 

682 0 FOR M=l TO 9 STEP 1 
20 683 0 A=0 

6840 FOR N=l TO 9 STEP 1 
6880 LET H=R (N) 
6890 H=-H 

6900 B=H*SIN(2*PI* (M-5) * (N-5) / 9) 
25 6910 A=A+B 

692 0 NEXT N 

693 0 A=A/9 

6935 LET MI (M) =A 

694 0 NEXT M 
30 6950 END SUB 

6960 SUB IFT (R() , I () , MR ( ) , MI () ) 
6970 DIM MRR ( 9 ) 
6980 DIM MR I (9) 
6990 DIM MIR (9) 
35 7000 DIM Mil (9) 

7010 FOR N=l TO 9 STEP 1 
7020 A=0 

703 0 FOR M=l TO 9 STEP 1 
7040 LET G=R(M) 
40 7050 B=G*COS (2*PI* (M-5) * (N-5) / 9) 
70 60 A=A+B 
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7070 NEXT M 

7080 LET MRR(N) =A 

7090 NEXT N 

7100 FOR N=l TO 9 STEP 1 
5 7110 A=0 

7120 FOR M=l TO 9 STEP 1 
713 0 LET G=R(M) 

7140 B=G*SIN(2*PI* (M-5) * (N-5) /9) 
715 0 A=A+B 
10 7160 NEXT M 

717 0 LET MR I (N) =A 
7180 NEXT N 

719 0 FOR N=l TO 9 STEP 1 
7200 A=0 
15 7210 FOR M=l TO 9 STEP 1 
722 0 LET G=I (M) 

7230 B=G*COS (2*PI* (M-5) * (N-5) / 9) 
724 0 A=A+B 
7250 NEXT M 
20 72 60 LET MIR (N) =A 
7270 NEXT N 

7280 FOR N=l TO 9 STEP 1 
7290 A=0 

7300 FOR M=l TO 9 STEP 1 
25 7310 LET G=I(M) 

7320 B=G*SIN(2*PI* (M-5) * (N-5) / 9) 

73 3 0 A=A+B 

7340 NEXT M 

7350 LET Mil (N) =A 
30 73 60 NEXT N 

7365 MAT MII=(-1)*MII 

73 75 MAT MR=MRR+MI I 

73 85 MAT MI=MIR+MRI 

7400 END SUB 
35 7410 SUB IFTZ (R ( } 7 I () , MR ( ) , MI () ) 

742 0 DIM MRR ( 9 ) 

743 0 DIM MRI (9) 
7440 DIM MIR (9) 
7450 DIM Mil (9) 

40 7460 FOR N=l TO 9 STEP 1 
7470 A=0 
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74 8 0 FOR M=l TO 9 STEP 1 

7490 LET G=R(M) 

7500 B=G*COS (2*PI* (M-5) * (-N+1) / 9) 

7510 A=A+B 

5 752 0 NEXT M 

753 0 LET MRR(N) =A 

7540 NEXT N 

7550 FOR N=l TO 9 STEP 1 

7560 A=0 

10 7570 FOR M=l TO 9 STEP 1 

7580 LET G=R(M) 

7590 B=G*SIN(2*PI* (M-5) * (-N+1) / 9) 

7600 A=A+B 

7610 NEXT M 

15 7620 LET MRI(N)=A 

7 63 0 NEXT N 

7640 FOR N=l TO 9 STEP 1 

7641 A=0 

7642 FOR M=l TO 9 STEP 1 
20 7650 LET G=I (M) 

7660 B=G*COS (2*PI* (M-5) * (-N+1) / 9) 

7 67 0 A=A+B 

7680 NEXT M 

7690 LET MIR (N) =A 

25 7700 NEXT N 

7710 FOR N=l TO 9 STEP 1 

7720 A=0 

773 0 FOR M=l TO 9 STEP 1 

7740 LET G=I (M) 

30 7750 B=G*SIN(2*PI* (M-5) * (-N+1) / 9) 

7760 A-A+B 

7770 NEXT M 

7780 LET Mil (N) =A 

7790 NEXT N 

35 7795 MAT MII=(-1)*MII 

78 00 MAT MR-MRR+MII 

7810 MAT MI=MRI+MIR 

7880 END SUB 

$ 

40 
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In an embodiment, the matrix inversion and/or the reiterative algorithms are used 
in combination with the Fourier Transform Algorithm. For example, an NMR image 
created by the Fourier Transform Algorithm may used as input for the first reiteration of 
the reiterative algorithm. 

5 

The Nyquist Theorem With The Determination Of The Spatial Resolution 

The derivation of Eq. (41) demonstrates that the system function behaves as a 
filter of the spectrum of the bulk magnetization function (NMR image). It is well known 

10 in the art of signal processing that if such a filter passes all frequencies for which an 
input function has significant energy, then the input function can be recovered 
completely from samples of the filtered function taken at the Nyquist rate. This premise 
embodies the Nyquist Sampling Theorem. The spectrum of the system function (Eq. (1)) 
is derived in APPENDIX V and shown in FIGURE lc. This function is a band-pass for 

15 all frequencies of the bulk magnetization function where k p and k z are comparable. 

Thus, the bulk magnetization function can be recovered by sampling the continuous RF 
field function given by Eq. (40) at the Nyquist rate, twice the highest frequency of the 
bulk magnetization function, in each spatial dimension over the sample space for which 
the function has appreciable energy. Sampling operations other than the present 

20 operation and the negligible error encountered by not sampling over the entire sample 
space are discussed in McC. Siebert [11] and the references therein disclosed which are 
all incorporated herein by reference. In the absence of noise, the spectrum of the bulk 
magnetization function can be completely recovered if the detector spacing frequency is 
equal to the Nyquist rate which is twice the highest frequency of the bulk magnetization 

25 function, and this represents the limit of resolution. However, the density of the detector 
spacing is limited by noise. The three-dimensional bulk magnetization map is a 
reconstruction from independent recordings at independent detector spatial locations 
relative to the voxels of the image space where two detector signals are independent if 
they are sufficiently spatially displaced from each other such that the difference in signal 

30 between the two detectors due to a test voxel is above the noise level. The resolution 
based on signal-to-noise arguments is discussed in the Contrast and Limits of Resolution 
Section. 
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Contrast and Limits of Resolution 



5 



The ability to visualize a structure in a noise-free environment depends, among 



other factors, on the local contrast C, which is defined as 




where / is the average background intensity and AI is the intensity variation in the 
region of interest. The main source of NMR image (also called magnetic resonance 
10 images (MRI)) contrast is T x and T 2 which depend on tissue types. The contrast may be 
increased by using RF pulse sequences to polarize the protons which encode T x 
(spin/lattice) and T 2 (spin/spin) relaxation time information in the data of the secondary 
RF field in a manner straightforward to those skilled in the NMR art. 



15 enhanced by, for example, subtracting part of the background or raising the intensity 
pattern to some power. Noise, however, represents a fundamental limitation on the 
ability to visualize structures. The signal-to-noise ratio, a basic measure of visualization, 
is determined by the ratio of the desired intensity variations to the random intensity 
variations, which are governed by the statistical properties of the system. The signal-to- 

20 noise ratio (SNR) is defined as 



where 0 7 is the standard deviation of the background intensity representing the rms 
value of the intensity fluctuations. The noise properties of the 4D-MRI imager involve 

additive noise only principally from thermal (y ) noise in the RF measurement circuits of 
25 the antennas of the detector array and from the fluctuations of the primary magnetic 



A feature of superconducting magnets is their extreme stability. The object to be 
imaged is magnetized with a highly stable magnet such as a superconducting magnet, hi 
this case, a magnetic field stability of 10~ 8 % over a month's time is feasible. Small 
30 antennas measure the RF signals as point samples without significant decrease in the 
signal to noise ratio relative to large antennas by using impedance matching while 



Contrast, is not a fundamental limit on visualization since it can be artificially 




(76) 



field. 
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minimizing resistive losses by using superconducting reactance elements, for example. 
In an embodiment, cross talk between antennas is ameliorated or eliminated by time 
multiplexing the signal detection over the array of antennas. External sources of RF 
noise can be ameliorated by placing the 4D-MRI scanner in a shielded room (Faraday 
5 cage). 

The quality of the image (i.e. the signal to noise ratio of the image) can be 
increased by repeating the reconstruction over multiple time points wherein each data set 
of a given time point is the set of matrices of the intensity variation over the sample 
space of the RF field of the bulk magnetization M of each voxel having the 

10 corresponding magnetic moment at that synchronized (common) time point. 

Sources of detector error are random, and the noise averages out as the number of 
detectors increases. Typically, the noise is suppressed by a factor of the inverse square 
root of the number of detectors. Thus, the effective limit of RF field detection is 
increased by a factor of the square root of the number of detectors for a constant SNR. 

15 The resolution depends on the extent that the field of the ring, shell, or sphere of 

dipoles differs from that of a single dipole at the center. The plot of the three cases of the 
field of a ring, shell, and a sphere of dipoles each of radius R and magnetic moment 
m = 10 4 Gem 3 given by Eq. (1.14), Eq. (11.17), and Eq. (IV. 16) of APPENDIX I, 
APPENDIX II, and APPENDIX IV, respectively, as a function of radius R where the 

20 position of the center of the ring, shell, and sphere relative to the detector is the point 

(0,0,10) is given in FIGURES 2, 4, and 6, respectively. The plot of the three cases of the 
field of a ring, shell, and sphere of dipoles of radius R = 0.2 cm and magnetic moment 
m = 10 4 Gem 3 given by Eq. (1.14), Eq. (11.17), and Eq. (IV.16) of APPENDIX I, 
APPENDIX II, and APPENDIX IV, respectively, as a function of the distance between 

25 the detector at the origin and the center of the ring, shell, and sphere at the points (0,0, 
z= 4 cm to z = 15 cm ) is shown in FIGURES 3, 5, and 7, respectively. From FIGURES 
2-7, it can be appreciated that the detector must be able to resolve three to four 
significant figures in order to reconstruct a map of 0.2 cm resolution with a field depth of 
15 cm. 

30 

Finite Detector Length 

The system function, h , of the Reconstruction Algorithm Section is the impulse 
response for a point detector. The following analysis will concern the impulse response 
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for an NMR detector which has finite dimensions. Consider a detector array 401 as 
shown in FIGURE 12 comprising multiple parallel planes 402 wherein each plane has a 
plurality of antennae coils 403 wherein each has an area A. The signal at any coil 403 is 
given by the integral of Eq. (40) over the area of the coil. This is the impulse response 
5 which is the system function which replaces h for a finite area detector. The Fourier 
Transform of this system function contains an argument of a product of the detector area 
A and the spatial frequency variables. Reconstruction could be performed as previously 
described in the Reconstruction Algorithm Section where this system function is 
substituted for the system function for a point detector. In the limit of zero area, the 

10 measurement is that of a point sample. Thus, another approach is to use the linearity of 
the superposition of RF signal to gain a higher signal to noise ratio advantage by taking 
the difference of relatively large signals as opposed to performing the measurements 
with miniature antennas. In this case, the set of n coils each of area A of each plane 
comprise a grid of n blocks each of area a«A formed by the absence of common area 

15 overlap of one or more of the coils. The signal of each element of the grid of area a can 
be solved by the set of linear difference equations of the signals of coils each of area A 
that corresponds to the noncommon area for the overlapping coils. The solution of the 
signals due to the grid elements can be obtain using a computer by matrix inversion. 
The resulting values represent the average signal for each grid center location. 

20 The effect of this data processing operation on the spectrum can be modeled as a sample 
and hold, where the voltages at the centers of the grid elements are sampled by 
multiplying by a picket fence of delta functions spaced s units apart which are convolved 
with d , a square wave function in the x-direction and the y-direction of width s units 
where the coordinates of Eq. (1) and FIGURE 9 are used. In the frequency domain, this 

25 data processing operation causes the spectrum of the signal function s to be multiplied 
by D, the Fourier Transform of the square wave function of width s units, to form 
function S * . If this multiplication does not multiply S , the Fourier Transform of the 
signal function, s , by zero for any frequency less than its bandwidth, then S can be 
recovered from S* by multiplying S* with the inverse of the Fourier transform of the 

30 sample and hold square wave function, a two dimensional sine function for the x and y- 
directions. This analysis applies to all axes in which direction the detectors have finite 
length. Furthermore, z-sampling is achieved by translating the array in the z-direction by 
interval distances at which points discrete signals are recorded or by using multiple 
parallel plane detector arrays spaced at the sampling interval along the z-axis. However, 
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if the signals are not sampled at discrete z-points, but each sample point is the integral 
resultant of the signal acquired continuously over a z-displacement of q units which is 
much greater than the dimension of the detector in the z-direction, then the 
corresponding sample and hold square wave has a width of q units. 

5 

It will thus be seen that the invention efficiently attains the objects set forth 
above, among those made apparent from the preceding description. Since certain 
changes may made in the above constructions without departing from the scope of the 
invention, it is intended that all matter contained in the above description or shown in the 
10 accompanying drawings be interpreted as illustrative and not in a limiting sense. In 

particular, other methods such as electronic and optical methods of detection of magnetic 
resonance are within the scope of the present invention. 

While the claimed invention has been described in detail and with reference to 
1 5 specific embodiments thereof, it will be apparent to one of ordinary skill in the art that 
various changes and modifications can be made to the claimed invention without 
departing from the spirit and scope thereof. 
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APPENDIX I 
Integration of a RING of Dipoles 

5 

Derivation of the Field Produced by a Ring of Magnetic Dipoles 

The z -component of the magnetic field due to a dipole or a loop of current of 
radius R with dipole moment m = inF? at the origin is given as follows: 

B -= 1(^ + /f J (L1) 

The z -component of the magnetic field at the position (x ? y,z) due to a ring of dipoles of 

m 

radius R with dipole density centered at the origin is 

2kR 



2 

B m 

z 2kR 



* (2 (z - z ' y - (x - x' ) 2 - - / ) 2 ><# 

i ((x-x0 2 + (^-y) 2 + (z-^) 2 ) 5/2 



wherein (x',y',z') is a variable which corresponds to the position of each individual dipole 
15 of the ring. The relationship between Cartesian coordinates and cylindrical coordinates 
with z 1 = 0 is 



r = ^jx 2 +y 2 +z 2 R = -Jx' 2 +y' 2 +z' 2 

x' = Rcos<l> (1.3) 
y = i?sin</> 

Substitution of Eq. (1.3) into Eq. (1.2) gives 

„, 2l i (lz 2 - (x - Rcos <j>¥-(y- sin^Y \<h 
B 2 B * W U (1 4) 

2 ^ o ((*-£cos0) +(j>-i?sin<£) +z 2 j 

20 Multiplying out terms gives 

m 2 p(2z 2 - (x 2 + R 2 cos 2 $ -2;ci?cos0)- (y 2 +i? 2 sin 2 -2.yj?sin4>))# 

B 2 = I — ~Z \5/2 Q-""0 

ln o ((x 2 + i? 2 cos 2 ^-2x/?cos^)+(^ 2 + i2 2 sui 2 0-2yi?sm0)+z 2 ) 

Substitution of Eq. (1.3) in the denominator and using the associative relationship and the 
trigonometric identity, cos 2 a + sin 2 a = 1 in the denominator gives 

m 2 *(lz 2 - (x 2 + R 2 cos 2 0 -2xRcos<£)- (v 2 +Rl sin2 0 -2>iRsin0))fy 
2 J (r 2 + i? - 2i?(x cos 0 + jysin 0)) 
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The denominator can be approximated using the following relationship 
(a + by = a*+nb(f- 1 

where 

a = r 4- R 



b = -2R(xcos<j) -t- vsin <f) 



(1.7) 

(1.8) 
(L9) 



and 



w = -5/2 (1.10) 
This case gives the far field where a «b and/or r»R. Using Eqs. (I.7-L10) gives 

(r 2 +i? 2 -2i?(xcos0 + y sin <j>jy 5/2 



+ ■ 



5i? 



J/2 



10 Substitution of Eq. (1. 11) into Eq. (1.6) gives 



(xcos</> +ysin$) 



B = 



T~ l 2z2 \ /2 1 2V« + /2 5 f 2 v/2 (* cos< ft +^sin0) <ty 

2?r o |_(r + *) ( r +i? ) J 



R 2 cos 2 <b 2xKcos<6 1,, 
+— t rr£r -— ttttt l# 



' 2 * JL(- 2 + * 2 f + <r 2 + *T (>- 2 + R T 1 

f 5i?x 3 cos (ft 5j?x 2 )^sin (/> 57? 3 xcos 3 <fr 
2tt J| ^ 5ig 3 < ycos 2 ^sin^ 10x 2 .fl 2 cos 2 0 lOxyi? 2 cos 0 sin ^> 



1 



L (r 2 + R 2 f 2 (? + R 2 f 2 (r 2 + R 2 f 2 j 



/77 

2n 



i 



y 



+■ 



R 2 sin 2 # 27igsin^> 



1 



r 2 + J R 2 ) 5 ' 2 ' (r 2 +i? 2 ) 5/2 (r 2 +if 2 

xcos^ 5Ry 3 sin (ft 5i? 3 xcos <ftsin 2 (ft 1 

-J(F7Ff' + P+* 2 ) 7/2+ (>- 2 +* 2 ) 7/2 1 

9 ^ J ' Si^ysin 3 ^ lOig 2 xycos^sin0 10R 2 y 2 sin 2 r i 

(r 2 + i? 2 ) 7/2 j 



+- 



L (r 2 +i? 2 ) ?/2 (r 2 +i? 2 ) 7/2 



(1.11) 



(1.12) 



The integration of Eq. (1.12) gives 

f Anz 2 2nx 2 tjR 2 + _ \0nx 2 R 2 1 



m l(r 2 + i? 2 Y /2 (r 2 + i? 2 ) 5/2 (r 2 + R 2 T (r 2 + R 2 f 2 \ 

b =— r y v 7 v ' v 7 j- an) 

2?r 2ttv 2 _ Trig 2 lOwfly , 

[ (r 2 + R 2 f 2 (r 2 + R 2 f 2 + (r 2 + R 2 f 2 j 
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Collecting terms gives 



5 

Ring of Dipoles: 

B z -m^ — - — ^572 + /2 V 2V /2 ^ ( L14 > 

10 

Eq. (1.14) demonstrates that the magnetic field of a ring of dipoles is not equal to 
that of a single dipole at the origin. The ring radius, R , appears in the denominator of 
each term. The first term is the magnetic field of a dipole at the origin only when the 

15 variable corresponding to the radius of the ring, R , is zero. The second and third terms 
are additional perturbations of the field of a dipole at the origin whose magnitude is a 
function of the radius of the ring. The second and third terms vanish only when the 
radius of the ring is zero. A ring of dipoles has a field that is cylindrically symmetrical. 
A shell, a cylinder, and a sphere of dipoles are the only other cases which have this 

20 symmetry. A cylinder is a linear combination of rings. Thus, the uniqueness of the 

dipole field is demonstrated by showing that it is different from that of a ring, a shell, and 
a sphere. The present result that the field of a dipole is different for that of a ring of 
dipoles as well as the same result in the cases of a shell and a sphere of dipoles shown in 
APPENDIX II and APPENDIX IV, respectively, demonstrate that dipole field is unique. 

25 All other fields are a linear combination of dipoles. Thus, the dipole is a basis element 
for the reconstruction of a NMR image. The resolution depends on the extent that the 
field of the ring of dipoles differs from that of a single dipole at the origin. The plot of 
the field of a ring of dipoles of radius R and magnetic moment m = 10 4 Gem given by 
Eq. (1.14) as a function of radius R is given in FIGURE 2. The position of the center of 

30 the ring relative to the detector is the point (0,0,10). The plot of the field of a ring of 
dipoles of radius R = 0.2 cm and magnetic moment m = 10 4 Gem 3 given by Eq. (1.14) 
as a function of the distance between the detector at the origin and the center of the ring 
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at the points (0,0, z= 4 cm toz = 15 cm) is shown in FIGURE 3. From FIGURES 2 and 
3 it can be appreciated that the detector must be able to resolve three to four significant 
figures in order to reconstruct a map of 0.2 cm resolution with a field depth of 15 cm. 
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APPENDIX II 

Integration of a Spherical Shell of Dipoles 

Derivation of the Field Produced by a Shell of Magnetic Dipoles 

The z -component of the magnetic field due to a dipole or a loop of current of 
radius R with dipole moment m = i%R 2 at the origin is given as follows: 



10 B =m 



I" 2z 2 -x 2 -y 2 "1 
[(x 2 +y 2 + z 2 f 2 _ 



(HI) 



The z -component of the magnetic field at the position (x ? y,z) due to a shell of dipoles of 

ffi 

radius R with dipole density r centered at the origin is 

4tvR 



2 

_ m 

z 4kR 2 



l ™ (2(z - 4 f - (x - xj -(y-y) 2 )R 2 sin $d$dQ 

I \ ((* -*') 2 + 0 - y ) 2 + (z- f) 2 ) 572 (IL2) 



wherein (x' ? y' ? z f ) is a variable which corresponds to the position of each individual dipole 
15 of the shell. The relationship between Cartesian coordinates and spherical coordinates is 

p=^ x 2 +y 2 + z 2 i?=v^ 2 +y 2 +^ 2 

x = p sin 0 cos 6) x } = Rsincj) cos 6 (IL3) 

j= p sin 0 sin 9 / = i?sin^sin0 

Substitution of Eq. (II. 3) into Eq. (II.2) gives 

271 



B = m 



'** fe(z-R cos <tf -(jc -flsin<ftcos0) 2 -(y- iteunftsine) 2 )?? 2 sm$d<t>d9 



47ri?2 oo ((jc-i?sin<ftcos0) 2 +(j-i?sin<ftsin0) 2 +(z-i?cos^) 2 j 

20 Multiplying out terms gives 

m 2 r"r (2z 2 -4zRcos<j> +2R 2 cos 2 <j> - x 2 +2xRsin<f>cos9 



(TI.4) 



B, = 



f f 

4%R 2 qo(x 2 - 2xR sin<ft cos 9 + R 2 sin 2 0cos 2 9 + y 2 - 2yRsm.$sm 9 



-R 2 sin 2 $ cos 2 9 -y 2 + 2j/i?sin0sin0 - R 2 sin 2 (ftsin 2 ^)^ 2 sin< 
+R 2 sin 2 0 sin 2 0 +z 2 - 2zi?cos0 + i? 2 cos 2 0) 5/2 

(E.5) 

Multiplying out the R 2 sin (j> term gives 
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m 2 } r (2.R 2 z 2 sin^ -4zi? 3 sin0cos(ft +2R 4 cos 2 (ft sin<ft - x 2 R 2 sin(ft +2xR 3 sin 2 (ft cos0 
z ~4nR 2 l\ (x 2 -2xRsia^cos9 + R 2 ^ ^cos 2 9+y 2 -2yRsm^sva9 

-R 4 sin 3 (ft cos 2 0 - y 2 R 2 sin (ft + 2yR 3 sin 2 <ft sin 9 - R 4 sin 3 (ft sin 2 0)f<ft d0 
+i? 2 sin 2 <ft sin 2 0 + z 2 - 2zR cos 0 + R 2 cos 2 tf> ) 5/2 

(H.6) 

Substitution of Eq. (II.3) in the denominator gives 

m 2 f p (2i? 2 z 2 sin (ft - 4zi? 3 sin 0 cos <ft + 2R 4 cos 2 <ft sin (ft - x 2 i? 2 sin (ft + 2xR 3 sin 2 (ft cos 9 
z ~ AtzR 2 \\ (p 2 -2xi2sin<ftcos9 + i? 2 sin 2 0cos 2 9-2^sin^sin0 

-i? 4 sin 3 (ft cos 2 0 - /i? 2 sin(ft + 2yi? 3 sin 2 (ft sin 9 -R 4 sin 3 (ft sin 2 9)d(j)d9 
+R 2 sin 2 ^ sin 2 9 - 2zR cos<f) + R 2 cos 2 (ft J 
5 (H.7) 

Using the associative relationship and using the trigonometric identity, 

cos 2 a + sin 2 a = 1 in the denominator gives 

m 2 ^ r (2.R 2 z 2 sin <ft - 4zR 3 sin <ftcos <ft + 2i? 4 cos 2 <f> sia<j> - x 2 R 2 sin<ft + 2xR 3 sin 2 (ft cos 6 
Bz ~ AkR 2 11 (p 2 +R 2 -2xRsm<j)cos9 

-R 4 sin 3 <ftcos 2 G-y 2 R 2 sin<j)+2yR 3 sin 2 4>sin0 -i? 4 sin 3 (ft sin 2 0)^0 
-2yi?sm<ftsin0 - 2zi?cos<ft) 5/2 

(H.8) 

10 The denominator can be approximated using the following relationship 

(a + bj = a n +nba"- 1 (II.9) 

where 

a = p 2 + i? 2 (11.10) 
b = -2xRsin <ft cos 9 — 2yR sin (ft sin 0 - 2zi?cos 0 (II. 1 1) 

15 and 

« = -5/2 ai-i2) 

This case gives the far field where a«b and/or p»R. Using Eqs. (II.8-II.12) gives 
(p 2 +R 2 —2xRsm<j) cos0 -2>>i?sin<ftsin0 -2zi?cos$) S/2 

i 5i? , , a . . . - ai-i3) 

s-; rryr + t 377 (;t sin (ft cos 0 + vsind>sin0 + zcos<ft) 

^ 2 +J R 2 / 2 (p 2 + i? 2 ) V 

Substitution of Eq. (HI 3) into Eq. (H.8) gives 
20 B z =— 3- I I — — —sTT + — — — ^7r(xsm(ftcos0+3;sin<ftsm0+zcos(ft) 
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4zR 3 sin$cos<£ 20zR 4 sin<f> cos<6 ^ . , n . . . _ ,\ 
_— 71 — „ 2 ns/2 71 — -fv7T^(*sin0cos0 +.ysin<£sin0 + zcos<£) 

(P +*) (p 2 + -K 2 ) 

2i? 4 cos 2 0sin(A 10.R 5 cos 2 0sin<i ✓ . . . . , . n ,x 

+ — ; \s/2 + — ; . 7 /2 (xsinffl cos0 + vsin<psm0 +zcosffl) 

(p 2 +R 2 f 2 (p 2 +R 2 f 

x 2 R 2 sm<b 5x 2 R 3 smS , . . „ . , . _ ,s 

-7-7 ^72—71 — 772 (xsin0cos0+^sin^sm0 + zcos^) 

(p 2 + R 2 ) (p 2 + R 2 ) 

IxR 3 sin 2 <£ cos 0 \0xR 4 sin 2 <£ cos0 / . . . . , . _ , 

+ — z -fj2 — — ~ T772 — (.tsinip cos0 +jsmip sm0 + zcos$ 

(p 2 + R 2 ) (p 2 +R 2 ) 



R 4 sin 3 <ftcos 2 0 5i? 5 sin 3 <ft cos 2 0 
(p 2 + R 2 f 2 ~ ( P 2 +R 2 f 2 



(xsin0cos0 + ,ysin0sin0 + zcos$) 



y R sm<b 5 y R sin0 / . , _ . , . n ,x 
— ■— r rhr — , -77T (xsmfflcosO + vsmfflsinfl + zcosa ) 

(p 2 +R 2 ) (p 2 + R 2 ) ■ 

2yR* sin 2 <b sin0 10vi? 4 sin 2 6 sin© / . , „ . , . n 

+-"7-; ; X 5/2 + — 7— rbn — (*sm0cos0 + >>sin0sin0 +zcos< 

(p 2 + R 2 ) (p 2 +R 2 ) 



R 4 sm 3 (bsm 2 6 5i? 5 sin 3 <£ sin 2 Q ( . , ■ . ■ n ±\ 

"77 — ? x5/2 - —77 — „-v7/i (* sin ^ cos 0 + jsin 0 sm0 + zcos<£) 
(p 2 + i? 2 ) (p 2 + irj 



1 



<fyrf0 (11.14) 



10 The integral of a sum is equal to the sum of the integrals. 



B 



m f 2 rfT 2ig 2 z 2 sin0 lORVsin^ . . . . , . . .J.. 
=7—3-1 I I \ , -• .,5/2 +77 — — w (xsm^cos0+j;sm0sm0+zcos0) ^0 

4?ri? [ooL(P (p +*) J 

2 r f[" 4zjR 3 sind>cos0 20zi2 4 sin#cos(/> / . , . . , . . ±^jxjn 
0 J 0 J |_ (p 2 + i? 2 ) (p 2 + R 2 ) J 

2 fV 2i? 4 cos 2 (isin^ 10i? 5 cos 2 d>sind> . . 



2n 



cos 0 +y sin0 sin 0 + z cos 

x 2 i? 2 sind> 5x 2 i? 3 sind> / . , - . , . _ ivLi Jrt 
Jn + ~o — "^-7^-(xsm0cos0 -f^smpsmO + ZCOS0J ^rfQ 



[p 2 + R 2 J' (p 2 +i? 2 ) 



2rew 



, , 2xi? 3 sin 2 d> cos 0 lOjci? 4 sin 2 6 cos0 ^ . . ^ . , . _ , N 
15 + I—77 — 3^72 — + — 77 — T-TtTI — (xsin0cos0+ysin4>sin0+zcos0; 

o J oL (P +R ) (P +R) 



dcjxiO 



2 r"^R 4 sin 3 (A cos 2 0 5R 5 sin 3 0 cos 2 0 , . , rt . , . . .J.j^n 
ll — 7 xs/2 — h — 7 .7/2 — (xsin^cos0 +3;sin0sin0 4-zcos^ aqxid 

o J oL (P +r 



■R 2 ) 5 ' 2 + (p 2 + R 2 ) 
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r r; y z i? 2 sin<f> 5y 2 R? sin<f> / . , . , . n • , 

J J 71 — ^^5/2 +7~; — ^7v772 (*sin$cos9 +^sinfsin0 + zcos#J |a$a0 



(p 2 +£ 2 )T ' (p 2 + R 2 ) 
2 ('{ t 2yR 3 sin 2 0 sin 9 10j^ 4 sin 2 0 sin 9 , . 



i 



2 f"J jg 4 sin 3 0 sin 2 9 5ig 5 sin 3 <ft sin 2 9 

" 0 J o J L (p 2 +^t + <p 2+ * 2 ? 



\7/2 



sin 0 cos 0 + j>sin <f> sin 0 + z cos 0) 

i i 

(xsin^cosG + >>sin0sin0 + zcos0)jrf06?0 \ 



(11.15) 



The integration of Eq. (11.15) gives 



B. =- 



to 



4^ 2 \(p 2 +R 2 f 2 3(p 2 +i? 2 ) ?/2 " 1 "3(p 2 +i? 2 ) 5/2 (p 2 +i? 2 ) 5/2 " 3(p 2 +i? 2 ) ?/2 



80?rz 2 i? 4 _ 87ri2 4 

.7/2 +" 



4roc 2 i2 2 _ 40ra; 2 i? 4 



4kR 4 



47vy 2 R 2 40 Try 2 R 4 

s.S/2 + 



4?ri? 4 



1 



3(p 2 + i? 2 ) 5 ' 2 3(p 2 + R 2 f 2 3(p 2 + i? 2 ) ?/2 3(p 2 + R 2 ) 5 ' 2 j 



(11.16) 



10 Collecting terms gives 



Shell of Dipoles: 



15 



B, = 



m 



r r 07 2 - r 
\akr?- z x 



AnR 2 [™ [ (, 



p 2 + 



„,v5/2 
i? 2 ) 



40 7c/2 4 


2z 2 -x 2 -/ 


3 





(11.17) 



20 Eq. (II. 17) demonstrates that the magnetic field of a shell of dipoles is not equal 

to that of a single dipole at the origin. The shell radius, R , appears in the denominator of 
the first term. The first term is the magnetic field of a dipole at the origin only when the 
variable corresponding to the radius of the shell, R 9 is zero. The second term is an 
additional perturbation of the field of a dipole at the origin whose magnitude is a 

25 function of the radius of the shell. The second term vanishes only when the radius of the 
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shell is zero. Thus, the dipole is a basis element for the reconstruction of a NMR image. 
The resolution depends on the extent that the field of the shell of dipoles differs from that 
of a single dipole at the origin. The plot of the field of a shell of dipoles of radius R and 
magnetic moment m = 10 4 Gem 3 given by Eq. (11.17) as a function of radius R is given 
5 in FIGURE 4. The position of the center of the shell relative to the detector is the point 
(0,0,10). The plot of the field of a shell of dipoles of radius R - 0.2 cm and magnetic 
moment m = 10 4 Gem 3 given by Eq. (11.17) as a function of the distance between the 
detector at the origin and the center of the shell at the points (0,0, z) is shown in FIGURE 
5. From FIGURES 4 and 5 it can be appreciated that the detector must be able to resolve 
10 three to four significant figures in order to reconstruct a map of 0.2 cm resolution with a 
field depth of 15 cm. 
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APPENDIX III 

Proof that the Field Produced by a Shell of Magnetic Dipoles is Different from that of a 
5 Single Dipole 



Consider Eq. (IL5). 

(2z 2 - 4zffcos <ft + 2R 2 cos 2 <j> - x 2 + 2xR sin<ft cos 0 



B = 



ff *f 

AnR 2 J J (x 2 - 2xi? sin<£ cos 9 + R 2 sin 2 0cos 2 6+y 2 - 2yRsin $ sin 0 

-i? 2 sin 2 0 cos 2 0 -j; 2 + 2.yi?sin 0sin0 - R 2 sin 2 0sin 2 e)/? 2 sin <^rf0 
+R 2 sin 2 $ sin 2 0 +z 2 - 2zRcos<f> + R 2 cos 2 fj' 2 

(in.i) 

10 The conditions for this integral to equal that of the field of a dipole at the origin are 

-4zRcos<t> + 2R 2 cos 2 <p + 2xRsin 0cos0 - R 2 sin 2 <£cos 2 0 + 2>i?sin^ sin0 - R 2 sin 2 ^sin 2 0=0 

(IIL2) 

and 

-2xi?sin (f) cos 0 + i? 2 sin cos 2 0 - 2j/#sin <j> sin0 + i? 2 sin 2 sin 2 0 - 2z£cos + J? 2 cos 2 = 0 
15 

(IIL3) 

Thus, Eq. (HL2) must equal Eq. (III.3). 
-4zRcos(j) + 2iJ 2 cos 2 $ + 2xRsin$cos0 - R 2 sin 2 <£cos 2 0 + 2yRsin<j> sin0 - R 2 sin 2 <£sin 2 0 
= -2xi?sin <j) cos 0 + i? 2 sin 2 0 cos 2 0 - 2j^i?sin <j) sin 0 + i? 2 sin 2 0 sin 2 0 - 2zi? cos <j) + R 2 cos 2 

20 (III.4) 
Therefore, 

-2zR cos <j> + i? 2 cos 2 <j>-2R 2 sin 2 0 cos 2 0 + 4x sin 0 cos 0 + 4yR sin 0 sin 0 - 2i? 2 sin 2 0 sin 2 0=0 

(IIL5) 

25 Using the associative relationship and the trigonometric identity, cos 2 a + sin 2 a = 1 as 
well as factoring out an R in Eq. (III. 5) gives 

-2zRcos (j> + R 2 cos 2 ^ - 2i? 2 sin 2 + 4x sin cos 0 + 4yR sin<£ sin0 = 0 (HI. 6) 
-2zcos <j) + 3R cos 2 - 2R + 4x sin<£ cos 0 + 4j>sin <j> sin 0 = 0 (IH.7) 

30 
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For x^0 ?t y^0, z*0, Eq. (HI. 7) is true only if R = 0 which proves that the field of a 
shell of dipoles centered on the origin is different from that of a single dipole at the 
origin. 

Eq. (III. 5) can not be integrated in closed form; however, the integral can be 
5 approximated so that the detector tolerances for a given image resolution can be 

determined. The denominator of Eq. (III. 5) can be approximated using the following 
relationship 

(a+by = a n +nba n ~ l (IIL8) 

where 

10 a = p 2 +R 2 (III.9) 

b = ~2xR sin 0 cos 0 - 2yR sin 0 sin 0 - 2zR cos (j> (III. 1 0) 

and 

77 = -5/2 (111.11) 
The result given by Eq. (II. 17) is 

15 B^T-jiM 2 — — 57T — — — -^f7T\} (IH.12) 



4kR 2 



Higher order terms of the approximation given by Eq. (III. 8) would contain mixed 
products of the coordinate variables which would increase the deviation of the derivation 
of the field a shell of dipoles from that of a single dipole at the origin. 
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APPENDIX IV 
Integration of a Sphere of Dipoles 

5 

Derivation of the Field Produced by a Sphere of Magnetic Dipoles 

The z -component of the magnetic field due to a sphere of dipoles is derived 
using the equation for a shell of dipoles with the substitution of the dipole density of the 

3 /72 

10 sphere ^ 3 where the radius of the sphere is R . The field of a shell of dipoles is given 

approximately by Eq. (II. 17) of the Derivation of the Field Produced by a Shell of 
Magnetic Dipoles Section 

m \ A J 2z 2 -x 2 -y 2 ^ 40nr^ 2z 2 -x 2 -y 2 \ mT ^ 

B z =— T {47ir 2 — ^TT" -— 71 } (IV.l) 

where the radius of the shell is r (The parameter r replaces the constant R of Eq. (11.17) 
15 as the radius of the shell.)- For small r , the first term of Eq. (IV.l) dominates and is 

used to approximately calculate the z -component of the magnetic field due to a sphere 

of z -oriented dipoles as follows: 

A sphere of dipoles is equivalent to the integral over concentric shells of dipoles 

each of radius r where 0<r< R. Thus, the 5 -component of the magnetic field due to a 
20 sphere of dipoles is the integral of the field of the shells given by Eq. (IV.l). 

B 2 = m 3 4n] ^f'^Zn ^dr (IV.2) 
2 4/3nR 3 J [_ (p 2 +r 2 f J 

Let 

- = tan0; dr = psec 2 QdQ (IV.3) 
P 

and using the trigonometric identity 
25 1 + tan 2 0 = sec 2 0 (IV.4) 

Eq. (IV.2) becomes 

It 1 

„ s 2 2 i\ m r An tan 0 2„™ ^ 

- rvR 3 0 P sec 9 
3 

Substitution in of 



WO 02/16956 



PCT/US01/25954 



96 



and 



tan0 = 



sec0 = 



sin 0 
cos 6 

1 



COS0 

into Eq. (IV. 5) gives 



A 

B 2 ^(lz 2 -x 2 -y 2 )-^— J -^sin 2 OcosGrfG 

— TlR 3 o P 

3 



Let 



u = sin 9; du = cos 0<i0 
Substitution of Eq. (IV.9) into Eq. (IV. 8) followed by integration gives 

m 4tvu 3 



B z =(2z 2 -x 2 -/>4 



3 

1 0 Substitution of Eq. (IV.9) into Eq. (IV. 1 0) gives 



From Eq. (IV.3) 



— = tanO; smQ = ■ 
P 



r 
P 



1 + 



f \ 2 
vp, 



Substitution of Eq. (IV. 12) into Eq. (IV. 11) gives 



15 



B 2 =(2z 2 -* 2 -/) 



2\ m 



Jfp 5 



Evaluation at the limits of the integral gives 



r 3 








1 + 




1) 


V 




J 



B z =(2z 2 -x 2 -y 2 ) 



2\ m 



i + i^ 



2 n 3/2 |o 
PJ ) 

2 \ 



Multiplication of Eq. (IV. 14) by 1 = • 



gives 



(IV.6) 
(IV.7) 

(IV.8) 

(W.9) 
(IV. 10) 



(IV.ll) 



(IV.12) 



(IV. 13) 



(IV.14) 
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B z =(2z 2 -, 2 -/) 



m 



R 1 



i + 



R 3 p 5 ' ' 



5/2 



(IV.15) 



Multiplying out the p 5 term in the denominator gives 



5 



Sphere of Dipoles: 




(IV.16) 



Eq. (IV.16) demonstrates that the magnetic field of a sphere of dipoles is not 
equal to that of a single dipole at the origin. The shell radius, R , appears in the 

15 denominator of the first term. The first term is the magnetic field of a dipole at the origin 
only when the variable corresponding to the radius of the shell, R , is zero. The second 
term is an additional perturbation of the field of a dipole at the origin whose magnitude is 
a function of the radius of the sphere. The second term vanishes only when the radius of 
the sphere is zero. Thus, the dipole is a basis element for the reconstruction of a NMR 

20 image. The resolution depends on the extent that the field of the sphere of dipoles differs 
from that of a single dipole at the origin. The plot of the field of a sphere of dipoles of 
radius R and magnetic moment m ~ 10 4 Gem given by Eq. (IV.16) as a function of 
radius R is given in FIGURE 6. The position of the center of the sphere relative to the 
detector is the point (0,0,10). The plot of the field of a sphere of dipoles of radius 

25 R = 0.2 cm and magnetic moment m = 10 4 Gent given by Eq. (IV.16) as a function of 
the distance between the detector at the origin and the center of the sphere at the points 
(0,0, z) is shown in FIGURE 7. From FIGURES 6 and 7 it can be appreciated that the 
detector must be able to resolve three to four significant figures in order to reconstruct a 
map of 0.2 cm resolution with a field depth of 15 cm. 
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APPENDIX V 



Fourier Transform of the System Function 



5 The system function, h(p 9 tf> 9 z) 9 in cylindrical coordinates is 

2z 2 — x 2 — y 2 2,z 2 — p 2 

h(p,(j),z)= 2 2 2 5/2 = 2 ^ 2,5/2 (V* 1 ) 

[x +y +z ] [p +z ] 
The spacetime Fourier transform in three dimensions in cylindrical coordinates, 
H(k p9 0 9 k z ) 9 is given [9] as follows: 

H(k p ,<S,k z ) = J J \h(p,<j>,z)exp{-i2x^k p pcos(<S> ~ ^) + K^>dpd^tbs (V.2) 

-oo 0 0 

10 With circular symmetry [9] 

OO 00 

H(k p ,k z )=27tj jh(p,z)J 0 (k p p)e-™pdpdz (V.3) 

0 -oo 

The Fourier transform of the system function is given by the substitution of Eq. (V.l) 
into Eq. (V.3). 

_ 2 2 

15 Consider the integral of Eq. (V.4) with respect to dp only. Factorization of h(p,(j) 9 z) 
gives 

-oof 2z 2 p p 3 1 

2« L I r 2 ^ 2 -,5/2 ~ r 2^ 2.5/2 P^p^P ( V ' 5 ) 
* LtP + ^ ] [p +Z ] J 

Consider the definite integral 

r t v+1 J v [at]dt a'z^K^JjazI 

Jo [^ + z 2 ]« +1 2«T[ii + 1] V ' } 

20 and the modified Bessel function of the third kind relationship, 

K„ v [x] = K u [x] (V.7) 

3 

The first factor of Eq. (V.5) is the same form as Eq. (V.6) with v = 0; u = — ? thus, 

2 % p 2z 2 (2rf'V 3 ' 2 [2 1/2 ]ttz 1/2 ^ /2 „ 

• 2z 2 (2n) j [p2 ^ z2]5/2 J 0 [k p p]dp= 2 a/2 r[5 P /2] ^- 3 / 2 [^ = r[5/2] P KsnW 

(V.8) 
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where K„ 3/2 [k p z] = K 3/2 [k p z] (Eq. (V.7)). The second factor of Eq. (V.5) can be made 

into the same form as Eq. (V.6) using the Bessel function of the first kind recurrence 
relationship 

5 Consider the second factor of the integral of Eq. (V.5). 

~ 27t f r £ *i W ( V - 10 > 
* [p +z ] 

Eq. (V.9) with v = 1 is 

^M+^M=-4M (V.H) 

J 0 W = -J 1 [x]-J 2 [x] (V.12) 

10 Let 

x=k p p (V.13) 

Substitution of Eq. (V.13) into Eq. (V.12) is 
2 

Substitution of Eq. (V.10) into Eq. (V.14) is 

3 r i 



•^^pP] = TT -At^pP] -J 2 [k p P] (V.14) 



15 



(V.15) 

The first factor of the right hand side of Eq. (V.15) is the same form as Eq. (V.6) with 

v = 1; u = — , thus, 
2 

20 - 2 < k p[ / +z 2 f 2 j ^ p]dp = K ~^ z] = - ' r[s/2] [ ^ 3 

(V.16) 

where iir_ 1/2 [A: z] = isTi /2 [A: p z] (Eq.(V.7)). The second factor of the right hand side of Eq. 

3 

(V. 1 5) is the same form as Eq. (V.6) with v = 2; w = — , thus, 
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[p 2 +z 2 f 2 " 2L " pl 



(2n)kl ,2 z U2 -nz in kT 
* K V2 [k p z] = 



2 3/2 r[5 / 2] " 1/2L ' V ^ J ~ p V2 y[5ii\ Kl,2 &] 



r 

5 \ 

j- 1 



(V.17) 

Combining the parts of the integration with respect to dp of Eq. (V.4) by adding Eq. 
(V.8), Eq. (V.16), andEq. (V.17) gives 

\2 v2 ~\v v2 m 



r[5 / 2] ' K3,l[ kpZ] r[5/2] Kin[kpZ] + ^^,^ K m tV] 



[2 1/2 ]T[5/2] 



The modified Bessel function of the third kind formulae is 



|- n -ji /2 « 



m!r[n + l -wj 



m=0 



Substitution of Eq. (V.13) into Eq. (V.19) with u = 1 is 



10 



K 3/2 [k p z] 



n 



2k p z 



1/2 



-,.r i i 
c 'L 1+ i^ r[3] J 



e 3Kz dz 
(V.18) 

(V.19) 
(V.20) 



Substitution of Eq. (V.13) into Eq. (V.19) with u = 0 is 



K V2 [k p z] = 



L2VJ 



-V 



Substitution of Eq. (V.20) and Eq. (V.21) into Eq. (V.18) is 

)7TZ 



fr^'v^r i i (2 i/2 )^ i7 v /2 ^ i/2 ^ /2 



J- j|_ H5/2] L 2 V J r £5/2] (2 1/2 )r[5/2]J 



15 



(V.21) 



r T /2 1 

(V.22) 



^ * -UK+Mz Z_Kl\3\ -UK+k p lz Z K -{jk^k^z , 7T , -[jk.+kjz 



-3/2 



■k e ^"""^ + — e 

r[5/2] p T[5/2]2 



T[5/2] 



r[5/2]2 9 



20 Collecting terms gives 

r[5/2] i J ' p 



rr[3] ~| _^ -uk+w 

i — 1 |Z 

L 2 J 



J^7T77^;i^[l + l/2] + ;^r :J -l |z _1 k "" " dz 



With T[3] =2 and r[5/2] = 3 /4tt 1/2 , Eq. (V.24) is 

3/2 

f — 4 [3 / 2] + [1 - 1JT 1 V^e^dz 



(V.23) 
(V.24) 

(V.25) 
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3/2 

\llf^ 3l2k ^ e ~ JKZdz (v - 26) 

2nk p [e~ k " z e- Jk ' 2 dz (V.27) 

4nk p £e~ k " z e- Jk * z dz (V.28) 

4Kk p [e- UWz dz (V.29) 
5 Integration of Eq. (V.29) with respect to dz gives 



4jlk J_^l_ e -t;w 



(V.30) 



p 

47T/J— ^— 1 (V.31) 



Multiplication of Eq. (V.31) by 

l^^f^l (V.32) 

10 gives 

A%K ^SK\ (V ' 33) 

The system function (Eq. (V.l)) is an even function; thus, the spacetime Fourier 
transform in three dimensions in cylindrical coordinates, H(k ,k z ), is given by taking the 



real part of Eq. (V.33) [8]. 



4nk 2 

15 mptk j = —Jt- (V.34) 

The spacetime Fourier transform in three dimensions in Cartesian coordinates, H{k p ,k 2 ), 
is 

2 , 7„ 2 



H[k x ,k y ,k z ]= ri 2 \ x i 2 (V.35) 



where the relationship between the wavenumbers and the spatial Cartesian coordinates is 
20 as follows: 
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lit 2% 
X~ z 



(V.38) 
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APPENDIX VI 



Derivation of S = HF® U(k 2 ) from Eq. (55). 



Ink 2 + ' i/2 +,i/2 r, , ,1 

S = HF® U[k x 1= S S Z Xn M e J ^ y ^ 

K P +K Z 2„ = 0 yn = -hl2 Xn^kll 

■ J-g) 47r/c p y V V y 

ilr It 2 4- lr 2 ^ ^ ^ A, h^i^h 

J^z Kp^ K zzn~0 y it = -h/2 x u =-hl2 

where the Fourier transform of u(z) = 1 for z > 0 and u(g) = 1 ybr z < 0 [8] is 

U(K) = ±8(k,)+± (VI.2) 

10 and where H is given by Eq. (52) and F is given by Eq. (49). The convolution integral 
of the second term of Eq. (VI. 1) is 

^ z„=o >.,=-/ 2 /2 *„=-/,/2 J _ K z |Ar p + \k z -K Z JJ 

Collecting exponential tenns gives 

15 (VI.4) 
Expansion of the denominator of the second term gives 

1 -/ 3 +/ 2 /2 +V 2 . r . , , tAttIc 2 



r 2 00 



(VI.5) 

20 The factorization of the denominator of the second term, k 2 z - 2k 2 K s + k 2 p + k 2 , using the 
quadratic formula is 



(VI.6) 
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2 

_ 2k z ±\Akl-Ak 2 p -Akl\ 
2 

k z = k z ±Jk p 
Substitution of Eq. (VI.6) into Eq. (VI.5) is 

1 -I, +U2 +IJ2 _ , 1 4 57 -t 2 " -, + ^' K * z » 

2 E X......^^'^-""^ f r _ t ' Tk _ k _ Jk fx, 

5 (VI. 7) 

The expansion of Eq. (VI. 7) by the method of partial fractions is 

-h + hl 2 +l i/ 2 



1 ~h +/ 2* z +t\t*> 



"f e+ >.». A+ £ + £ "1 

-i 13 K z-K+J k P K z -k z -jk p ] 



(VI.8) 



7 

The factors A, B, C of Eq. (VI.8) are determined as follows: 

a\k 2 z -2k z K z + k 2 p + k\ ]+ B[k z ][k z - k z -jk p ]+ C[k z Jk z -k z + jk p ]= 1 (VI.9) 
10 Let k z = 0. Then, Eq. (VI.9) gives 

4*J+#]=1 (VI.10) 

Thus, 

^^F] (VL11) 

Let k 2 = £ 2 + jk p . Then Eq. (VI.9) is 

^ Pi + JK T ~ + A 1+ ^ + \ % +JK lK +jk p - k z -jk p ] 

jgltJLJ J L A J (VI.12) 

+ C]k z + jk p Jk z + jk p —k z + jk p ]= 1 
The term in B vanishes. Substitution of Eq. (VI. 11) into Eq. (VI.12) gives 

^rjFp +2JKkp ~ kp ~ 2k * ~ 2jKkp +kt <> +^] +c [ 2 -M> - 2k l]= 1 (v 113 ) 

where the term in A also vanishes. Solving for C gives 

C - [-.2^+W/| <™> 

20 Let K 2 = k 2 ~jk . Then, Eq. (VI.9) is 
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2 co 



10 Thus, 



15 



is 



[-2% - 2jk 2 k p ]_{ j[-2k 2 ] 



(VI.15) 



A Ik -jk] - 2K - jk ]+ k + k \ b[k - jk p jk 2 - jk p -k,- jk p ] 
+ c% - jk p \k 2 ~jk p -k 2 + jk p y i 

The term in C vanishes. Substitution of Eq. (VI. 1 1) into Eq. (VI.15) gives 
[k'lk 2 ^ ~ 2Jk * k r ' e o~ 2 % +2Jk * k » + % +%]+B[-2jk 2 k p -2k 2 p ]= 1 (VI.16) 
where the term in A also vanishes. Solving for B gives 

B = T =, (VI.17) 

\-2kl -2jk 2 k p ] 

Substitution of A , B , and C into the convolution integral (Eq. (VI.8) gives: 

Z 2 „=o y, t =-h/2 x n =-hll 



_ ^_ I (VL18) 

I J i[<,-K+Jk p J r 2k p 2 -2jk 2 k p J ^ 



The first convolution integral is of the form given by Mc Seibert [8]. 
2 

— <=> sgnz n ; sgnz n = -1 where z n < 0 (VI. 19) 



4tt& 2 °° r e +jKzZ " 2nk 2 
A change of variable in the second convolution integral given by letting 

k, = -JK > dK * = ~ dk * ( VL21 > 

in 

Ank 2 % e^' 2 " 



— 2 — J —r^rrl- dk A where z « <0 ( VL23 > 
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| CVL24) 



e -^ |Z »' where z„ <0 (VI.25) 



+A* P ] 

A change of variable in the third convolution integral given by letting 
5 K z = -k z +jk p . dK z = -dk z (VI.26) 



m 



\-2kj +2Jk 2 k p ]i j[ K , -k, -jk p J ' 
Ink 1 % 2 

r : e =, f ££ -F T~\-dkA where z n <0 (VI.28) 



is 



10 



1 f — dk z (VI.29) 

[-2^+2^^ p ]i 2A 



ie + * pkl where z„ <0 (VI.30) 



2f p 2 -A* P ] 

Combining Eqs. (VI.20), (VI.25), and (VI.30) gives the convolution of Eq. (VI. 1). 



-'3 +/ 2' 2 +'i' 2 r if Trfr 2 1 1 irk 2 . 1 il 



(VL31) 
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APPENDIX VII 

Derivation of the Inverse Transform of Eq. (68) to Give Inverse Transform 1, Eq. (69). 

5 

l*z„=0 y n ~-h\2 x n =-k\2 L *p L J JJ 

(VIL1) 

C oo co ^ 

I ~ f f -k p \z„\ -r Pr Ht 7T ilr r. tt I 



1A S S 

10 o 



*n *» 



2tt J J ^J'J/.^p^d^d^ 

-oo 0 
-oo 0 

CO CO 



-oo 0 

CO CO 



I -72* J \e k ^J 0 ]k p p^k p k z e^dk z 

I -00 0 



where z„ < 0 (VII.2) 

Consider the definite integral: 

r -«* , r, i v+i j 2a[2feJ , rJv + 3/2] n _ m ,. 
J e "'jXbtY dt = \ J ' ^ a > 0 (VII.3) 
; [a + & 2 J 7T 1 ' 2 

In the case of the first integral of Eq. (VII.2), the parameters of Eq. (VH3) are 
v = 0 
t = k p 

15 , (VII.4) 

b = p 

a = \ Z n\ 

Thus, 

oo 
0 

where 

2 W 3 / 2 ] y_ (vii 6) 



o |5« + p J 



WO 02/16956 PCT/US01/25954 

108 

Consider the following relationship of the modified Bessel function of the third kind: 



+y*[v+i] r -y-[v+i] r -j- 

nJ v \z\ = e 2 K v ^ze 2 j + e 2 .K^ze 2 



o 

= — J |[cosfc p z„]+ jsin[/fc p z„]f cos[fc p z J-jsinf^ J^[/c p p)//c p (Vn.12) 



-2p" 2 



1 5 The transform is given by B ateman [12]. 
Thus, 



(vn.7) 



In the case that v = 0 , Eq. (VTL7) is 

tvJ 0 [£ p p]= jK 0 [jk p p]- jK a [-jk p p] (VII. 8) 

5 Thus, 

j" e +k » Kl J 0 \k pP Jfc p = J - \jK 0 [jk p p]-jK 0 [-jk pP $k p (VH.9) 

o o n 

= \ J'^ir^ l/VK -4 1 e^\[-jk pP yk p 

0 JL ' t 0 

(VII. 10) 

With a change of variable as follows: 
10 First 
Integral: Second Integral: 

k p =jk p k p =-Jk p 

dk p =jdk p dk p = -jdk p 

and an expansion of the complex exponential factor, Eq. (VE.10) becomes 



(VH.11) 



(VE.13) 



) e k ^J a ]k pP yk p = * (Vn. 14) 

The Fourier transform relationship between a function /[x] and its Fourier transform 
20 g-|>] given by McC. Siebert [8] is 
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xf[*]^-~gm (VH.15) 
J dk 

Consider the following Inverse Fourier transform integral from Eq. (VH2) 

J k/^j 0 [k pP ]ik p (vn.16) 

0 

Eq. (VII. 16) is equivalent to the following Inverse Fourier transform 

5 V^&fep]} ( VII - 17 > 

FromEq. (VH.15) with 



Eq. (VH.16) is 



(vn.i9) 



[4 

1 0 Taking the imaginary part because Eq. (VII. 1 9) has odd symmetry gives 

J k/^j 0 [k pP p p = r2 ~^"L 2 ; z n <o (vn.20) 

o |f » + P J 

In the case of the third and fourth terms of Eq. (VIL2), consider the following Fourier 
transform given by Gray [13]: 

J e"J a [bx]dx = 1 ^ (VII.21) 

15 Thus, 

J e-^j e [k pP }k p = 1 2f2 ; z„ > 0 (VH22) 

The Inverse Fourier transform with respect to k z is given by McC. Siebert [8]. In the case 
of the first two terms of Eq. (VIL2), 

oo 

\ e +Jk ' z dk z = 8[z] (VE.23) 

— oo 

20 In the case of the third and fourth terms of Eq. (VII.2), 

oo 

\jk z e +Jk ' z dk 2 =8[z'\ (VE.24) 
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Combining transforms from Eqs. (VII.6), (VTL20), and (VII.23) in the case of the first 
and second terms of Eq. (VIL2) and combining transforms from Eqs. (VII. 14) (VII.22) 
and (VII.24) in the case of the third and fourth terms of Eq. (VII.2) gives: 



F 



l°*„=o Vn—hn ^=-/i/2 L k p L -tlj 



5 

where z n <0 

(VH.25) 

Substitution of Eq. (48) into Eq. (VII.25) gives 

r»,-<> ^=-'2/2 ^=-4/2 L ^p 1 - 4IJ 

= "2 £ S X m ,, e , m S[x- ni k,y-n 2 k,z-n 3 k] 

5 « 3 =0 « 2 = -/ 2 /2i n^-ljlk 

10 

where z n <0 

(VH.26) 

Taking the derivative given by the doublet function gives 

[O^, = o y tt =-hl2 x„=-U2 L L J JJ 

+ +/ 2 /2& + ^2* f "] 

^7(3 = 0 n 2 =-hl2k ny=-kl2k L If" ^ J J 

15 

Inverse Transform 1 (VTL27) 
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CLAIMS: 

1 . Apparatus for providing a multidimensional image of inanimate or animate objects 
comprising: 

5 a magnetization source for producing a magnetic field to magnetize a volume of 

an object to be imaged, the volume comprising a plurality of voxels; 

a radiation source for applying a radiation field to the object to be imaged when 
the object is disposed in the magnetic field to produce a secondary radiation field 
emanating from the object; 
10 an output signal detector for producing output signals in response to said 

secondary radiation at a plurality of spatial locations outside of the object as a 
function of time; 

a first processor for determining a plurality of Fourier components, each having 
the same frequency, an intensity and a phase angle; 
15 a second processor for associating the Fourier components due to each voxel of 

the object to be imaged by phase to form a set corresponding to a spatial variation of 
intensity of the secondary radiation due to each voxel at said plurality of spatial 
locations outside of the object; 

a third processor for converting each said set of components into a voxel 
20 location; and 

an image processor for producing an image based on the secondary radiation 
intensity from each voxel and the location of each voxel. 

2. The apparatus according to claim 1, further comprising a display for displaying the 
25 image. 

3. The apparatus according to claim 2, characterized in that said display comprises at 
least one selected from the group consisting of printers, cathode ray tube displays, 
liquid crystal displays, plasma screens, three dimensional modelers, laser monitors, 

30 projections monitors and holographic displays. 

4. The apparatus according to any one of the preceding claims, characterized in that one 
processor handles the processing functions of at least two of said first processor, 
second processor, third processor and image processor. 
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5. The apparatus according to any one of the preceding claims, characterized in that one 
processor handles the processing functions of at least three of said first processor, 
second processor, third processor and image processor. 

5 6. The apparatus according to any one of the preceding claims, characterized in that the 
third processor comprises a reconstruction processor. 

7. The apparatus according to any one of the preceding claims, characterized in that the 
image is a magnetic resonance image. 

10 

8. The apparatus according to any one of the preceding claims, characterized in that the 
magnetization source comprises a superconducting magnet. 

9. The apparatus according to any one of the preceding claims, characterized in that the 
15 radiation source comprises a radio frequency generator and transmitter. 

10. The apparatus according to claim 9, characterized in that the radio frequency 
generator and transmitter comprises an antennae to excite protons in a magnetized 
volume. 

20 

11. The apparatus according to claim 10, characterized in that the antenna comprises two 
saddle coils driven in parallel to provide a rotating radio frequency field. 

12. The apparatus according to any one of the preceding claims, characterized in that the 
25 source of radiation produces a component of secondary radiation which is transverse 

to a direction of magnetization of a volume. 

13. The apparatus according to claim 12, characterized in that the secondary radiation 
which is transverse to the direction of the magnetization of the volume comprises 

30 radio frequency dipole radiation. 

14. The apparatus according to claim 13, characterized in that the radio frequency dipole 
radiation comprises near field radiation. 
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15. The apparatus according to any one of the preceding claims, characterized in that the 
output signal detector comprises a plurality of detectors. 

16. The apparatus according to claim 15, characterized in that each of the detectors 
5 comprises an associated antenna. 

17. The apparatus according to claim 16, characterized in that the plurality of detectors 
comprises a three dimensional array of antennas. 

10 18. The apparatus according to claim 17, characterized in that each antenna of the three 
dimensional array of antennas responds to a secondary radio frequency field as a 
function of time. 

19. The apparatus according to claim 18, characterized in that the plurality of detectors 
15 are located at spatial locations outside of an object to be imaged. 

20. The apparatus according to any one of the preceding claims, characterized in that the 
output signal detector produces output signals in response to secondary radiation at a 
plurality of spatial locations as a function of time samples at a Nyquist rate in time 

20 over the proton free induction decay. 

21. The apparatus according to any one of the preceding claims, characterized in that the 
first processor further comprises an analog to digital converter to digitize signals of 
secondary radiation. 

25 

22. The apparatus according to any one of the preceding claims, characterized in that the 
first processor comprises a time Fourier transform processor to convert secondary 
signals from the signal output detector over time into Fourier components, each 
having the same frequency, an intensity and a phase angle. 

30 

23. The apparatus according to any one of the preceding claims, characterized in that the 
plurality of Fourier components each having the same frequency, an intensity and a 
phase angle are such that the phase of each component is unique. 
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24. The apparatus according to claim 23, characterized in that the unique phase of each 
component is provided by a combination of the angle 6 suspended between each 
spatial location at which the secondary field is detected and a radial vector, a vector 
from each voxel to each spatial location, and the angle $ due to a separation distance 

5 r between each voxel and each spatial location given by a wavenumber of detected 

field k times r . 

25. The apparatus according to any one of the preceding claims, characterized in that 
said first processor uses a first component having a phase angle, calculates the phase 

10 angle as a function of spatial position of a first spatial position relative to any other 
spatial location and identifies a component at each spatial location having the 
calculated phase angle. 

26. The apparatus according to claim 25, characterized in that a separate detector 
15 samples each separate spatial location. 

27. The apparatus according to claim 26, characterized in that each detector comprises a 
point detector. 

20 28. The apparatus according to any one of the preceding claims, characterized in that the 
source of radiation comprises at least one radio frequency generator selected from the 
group consisting of klystrons, backward wave oscillators, Gunn diodes, and traveling 
wave tube amplifiers. 

25 29. The apparatus according to any one of the preceding claims, characterized in that the 
third processor uses a spatial variation of the secondary radiation over a sample 
volume comprising a plurality of spatial locations at which the radiation is detected 
to determine a location of the voxels. 

30 30. The apparatus according to claim 29, characterized in that the location of each voxel 
is determined independently by the third processor. 
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3 1 . The apparatus according to claim 2, wherein the image processor and display can 
rotate and display an image from a plurality of desired perspectives. 

32. The apparatus according to claim 2, characterized in that the image processor and 
5 display can display a reduced dimensional image. 

33. The apparatus according to any one of the preceding claims, characterized in that the 
apparatus provides a two dimensional or three dimensional nuclear magnetic 
resonance image of inanimate or animate objects. 

10 

34. The apparatus according to claim 33, characterized in that the output signal detector 
comprises a nuclear magnetic resonance detector for detecting a time dependent 
transverse radio frequency field component as a function of time of over a three 
dimensional sample volume outside of an object. 

15 

35. The apparatus according to claim 34, characterized in that secondary radiation 
emanating from a magnetized object comprises the transverse, time harmonic nuclear 
magnetic resonance magnetic dipole radiation. 

20 36. The apparatus according to claim 34, characterized in that the time dependent 

transverse radio frequency field component comprises a harmonic oscillation of each 
radio frequency dipole that is equivalent to a dipole rotating in a transverse plane. 

37. The apparatus according to claim 34, characterized in that the apparatus can generate 
25 an image from a three-dimensional map of a transverse resonant radio frequency 
magnetic flux external to an object and said nuclear magnetic resonance detector 
comprises a plurality of detectors at spatial locations, and wherein 

a nuclear magnetic resonance signal of each voxel at a given detector gives rise to 
a Fourier component with a unique phase angle relative to a Fourier 
30 component of other voxels of the object at that detector; 

a set of Fourier components that correspond to the nuclear magnetic 
resonance signal of a given voxel over the detectors is determined; and 

an intensity variation of a transverse radio frequency field over space is 
used to determine a coordinate location of each voxel. 
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38. The apparatus according to claim 37, characterized in that the each detector is 
selectively responsive to a component of the secondary field. 

5 39. The apparatus according to claim 38, characterized in that a transverse radio 

frequency field is a near field which is a dipole that serves as a basis element to form 
a unique reconstruction. 

40. The apparatus according to claim 39, characterized in that the each detector is 
10 responsive to a component of radio frequency in a transverse plane. 

41. The apparatus according to claim 40, characterized in that each detector is 
responsive to a component of the radio frequency dipole field parallel to the radio 
frequency dipole. 

15 

42. The apparatus according to claim 41, characterized in that each detector is 
responsive to the component of the radio frequency dipole field parallel to the radio 
frequency dipole as a function of time and detects a maximum intensity of the signal 
during each time harmonic cycle. 

20 

43. The apparatus according to any one of the preceding claims, characterized in that the 
third processor determines a location of a voxel from signals recorded of a transverse 
dipole radio frequency field. 

25 44. The apparatus according to claim 43, characterized in that a geometric system 

function corresponding to a dipole determines a spatial intensity variations of the 
radio frequency field. 

45. The apparatus according to claim 44, characterized in that the system function is a 
30 band-pass when k p = £ z? the frequency of spatial variation of the recorded secondary 

radiation in the rho direction is equal or comparable to that in the z direction of the 
three dimensional space comprising said plurality of spatial locations. 
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46. The apparatus according to claim 45, characterized in that each voxel is 
reconstructed independently in parallel with other voxels using each unique set of 
components associated by phase such that the scan time is no greater than the nuclear 
free induction decay time. 

5 

47. The apparatus according to claim 46, characterized in that all voxels comprising an 
image are reconstructed in parallel. 

48. The apparatus according to claim 46, characterized in that the apparatus can generate 
10 high resolution, three-dimensional, real-time anatomical images and images based on 

physiological parameters with little or no deterioration from motion artifact. 

49. The apparatus according to claim 41, characterized in that at each point in time, each 
radio frequency dipole rotating in the transverse plane is directed at an angle 6 

1 5 relative to the direction of detection of each detector. 

50. The apparatus according to claim 49, characterized in that the phase angle 6 of the 
radio frequency dipole relative to the direction of detection axis of the detector gives 
rise to a phase angle term e~ l& . 

20 

51. The apparatus according to claim 50, characterized in that at each point in time, each 
radio frequency dipole rotating in the transverse plane with a separation distance r 
between each voxel and each spatial location at a detector gives rise to an angle <j> 
given by a wavenumber of the detected secondary radiation k times r directed at the 

25 angle 9 relative to a direction of detection of each detector. 

52. The apparatus according to claim 51, characterized in that the phase angle (j> of the 
radio frequency dipole relative to a direction of detection axis of a detector gives rise 
to a phase angle term directed at angle 6 relative to the direction of detection of 

30 each detector and wherein the angle 6 gives rise to a phase angle term e~ lQ . 

53. The apparatus according to claim 52, characterized in that a sum of the phase angles, 
kr and 0 , is unique for each voxel at each detector. 
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54. The apparatus according to claim 37, characterized in that the unique phase of each 
component is provided by a combination of 

angle 0 suspended between each spatial location at which a secondary nuclear 
5 magnetic resonance field is detected and a radial vector from each voxel to each 

spatial location, and 

angle $ due to a separation distance r between each voxel and each spatial 
location given by the wavenumber of detected nuclear magnetic resonance field k 
times r . 

10 

55. The apparatus according to claim 37, further comprising a processor, characterized 
in that a set of Fourier components that correspond to a nuclear magnetic resonance 
signal of a given voxel over detectors is determined by the processor using a position 
of a first detector relative to a different detector to calculated the phase angle of the 

15 different detector relative to the first detector. 

56. The apparatus according to claim 55, characterized in that the processor determines 
the set of Fourier components that correspond to the nuclear magnetic resonance 
signal of a given voxel over the detectors by using the position of each first detector 

20 relative to a different detector to calculated the phase angle of the different detector 
relative to the first detectors, and the process is repeated over all of the detectors to 
give the set of intensities of the nuclear magnetic resonance signal over the plurality 
of detectors due to each voxel. 

25 57. The apparatus according to claim 56, characterized in that the set of intensities of the 
nuclear magnetic resonance signal over the plurality of detectors due to each voxel is 
processed to determine a coordinate location of each voxel. 

58. The apparatus according to any one of the preceding claims, further comprises coils 
30 for providing 

B Q , a steady primary field, 
G x , a field gradient in an x-axis direction, 
G y , a field gradient in a y-axis direction, 
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a radio frequency field, and 

G z , a field gradient in a z-axis direction. 



59. The apparatus according to claim 58, associated amplifiers for driving the coils and 
5 associated control circuits for controlling the coils which provide T x and T 2 pulse 

sequences to alter an intensity of the secondary radio frequency field as a function of 
time to provide corresponding signals to be processed. 

60. The apparatus according to claim 9, characterized in that nuclei present in object are 
10 subjected to an additional rotating radio frequency field Hj , which is synchronous 

with precession and their magnetic moments, and magnetization M is caused to 
precess about H t and rotate away from the primary field H 0 by an angle ^ in a 
coordinate frame which rotates at the Larmor frequency, the precession about H x 
continues as long as H 1 exists, and a final value of which depends on the strength 
15 of H 19 which determines the precession rate, and the time for which it is turned on is 
such that 0 = 90° that a dipole is perpendicular to a detector array that is responsive 
to the nuclear magnetic resonance radiation at a plurality of spatial locations external 
to the object. 

20 61. The apparatus according to any one of the preceding claims, characterized in that the 
magnetizing source is constructed to provide a magnetic field which permeates the 
object to be imaged and is confined to a volume to be imaged, and the confined 
magnetic field limits the source of signal to the volume of interest, characterized in 
that the volume to be reconstructed is limited to the magnetized volume which sets a 

25 limit to required computations and reduces end effects of signal originating outside 
of edges of the output signal detector. 

62. The apparatus according to claim 61, characterized in that the apparatus is 

constructed to provide a magnetizing field which limits an imaged volume of an 
30 object by providing a range of Larmor frequencies and data is rejected at a threshold 
Larmor frequency. 
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63. The apparatus according to any one of the preceding claims, characterized in that the 
apparatus is constructed to provide a magnetizing field which limits an imaged 
volume of an object by providing a range of Larmor frequencies and data is rejected 
at a threshold Larmor frequency. 

5 

64. The apparatus according to any one of the preceding claims, characterized in that the 
output signal detector comprises detector elements disposed on a first axis and 
having corresponding output signals. 

10 65. The apparatus according to claim 64, characterized in that the output signal detector 
includes detector elements disposed on a second axis and having a corresponding 
output signal, and the detectors on the second axis are responsive to the 
corresponding oriented field components of the secondary radiation. 

15 66. The apparatus according to any one of the preceding claims, characterized in that 
the third processor comprises at least one reconstruction processor selected from the 
group consisting of a detector processor for processing the detector output signal 
with an inverse matrix, a reiterative processor, a time Fourier transform processor, a 
sine function filter processor for detecting the detector output signals, a space and 

20 time Fourier transform processor, and a compensating processor for compensating 
for a distance from the detector along a detector responsive axis. 

67. The apparatus according to any one of the preceding claims, characterized in that the 
signal output detector further comprises a line or plane of detectors and a drive 

25 mechanism that moves the detectors to sample an external secondary radiation over a 
three-dimensional volume. 

68. The apparatus according to claim 67, characterized in that the secondary radiation is 
a radio frequency magnetic field due to nuclear magnetic resonance of nuclei in a 

30 magnetic field produced by the magnetization source and excited by a radiation field 
from the source of radiation; 



69. 



The apparatus according to claim 68, characterized in that the nuclear magnetic 
resonance active nuclei are protons. 
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70. The apparatus according to claim 19, characterized in that each antenna of the three 
dimensional array of antennas are impedance matched. 

5 71 . The apparatus according to claim 70, characterized in that each antenna of the three 
dimensional array of antennas further comprises superconducting reactance elements 
which minimize resistive losses. 

72. The apparatus according to claim 19, characterized in that cross talk between 
antennas is ameliorated or eliminated by time multiplexing the signal detection over 
the array of antennas. 

73. The apparatus according to claim 19, characterized in that the signal output detector 
comprises a synchronous detector wherein the radio frequency field is sampled 
synchronously so that a phase at any given detector may be related to that at any 
other detector. 

74. The apparatus according to claim 19, characterized in that the signal output detector 
comprises a synchronous detector wherein the radio frequency field is sampled at 
known times so that a phase at any given detector may be related to that at any other 
detector. 

75. The apparatus according to claim 19, characterized in that each antenna of the three 
dimensional array of antennas responds selectively to a component of the secondary 
radio frequency field at a spatial location outside of the body as a function of time. 

76. The apparatus according to claim 19, characterized in that each antenna of the three 
dimensional array of antennas responds selectively to a parallel component of the 
secondary radio frequency field at a spatial location outside of an object as a function 
of time. 

77. The apparatus according to claim 19, characterized in that each antenna of the three 
dimensional array of antennas responds selectively to a perpendicular component of 
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the secondary radio frequency field at a spatial location outside of an object as a 
function of time. 

78. Apparatus for providing a multidimensional image of inanimate or animate objects 
comprising: 

means for producing a magnetic field to magnetize a volume of an object to be 
imaged, the volume comprising a plurality of voxels; 

means for applying a radiation field to the object to be imaged when the object is 
disposed in the magnetic field to produce a secondary radiation field emanating from 
the object; 

means for producing output signals in response to said secondary radiation at a 
plurality of spatial locations outside of the object as a function of time; 

first processor means for determining a plurality of Fourier components, each 
having the same frequency, an intensity and a phase angle; 

second processor means for associating the Fourier components due to each voxel 
of the object to be imaged by phase to form a set corresponding to a spatial variation 
of intensity of the secondary radiation due to each voxel at said plurality of spatial 
locations outside of the object; 

third processor means for converting each said set of components into a voxel 
location; and 

image processor means for producing an image based on the secondary radiation 
intensity from each voxel and the location of each voxel. 

79. An apparatus for providing a two dimensional or three-dimensional nuclear 
25 magnetic resonance image of an object comprising: 

means for generating an image from a three-dimensional map of a transverse 
resonant radio frequency magnetic flux external to an object having a volume to be 
imaged, the volume comprising a plurality of voxels; and 

processor means for determining the location of each voxel which is a source of 
30 the radio frequency field through spatial variations of an intensity of the nuclear 

magnetic resonance field of a set of components associated by phase. 

80. An apparatus for providing a two dimensional or three-dimensional nuclear 
magnetic resonance image of an object comprising: 



10 



15 
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a magnetizing source for magnetizing a volume an object in a first direction, the 
volume comprising a plurality of voxels; 

means for causing nuclear magnetic resonance magnetization to rotate into a 
plane transverse to the first direction; 

an output signal detector for detecting the nuclear magnetic resonance field 
rotated into a transverse plane external to the object at a plurality of spatial 
locations; 

a first processor for determining a set of components of spatial variation of the 
nuclear magnetic resonance field over a sample volume by association of the 
components by phase; 

. a second processor for determining a location of each voxel which is a source of 
the radio frequency field through each set of components representative of spatial 
variations of an intensity of the nuclear magnetic resonance field over the sample 
volume, and 

an image processor for displaying the voxel locations representative of 
anatomical images of the object based on selected physiological 

81 . An apparatus for providing a two dimensional or three-dimensional nuclear 
magnetic resonance image of an object comprising: 

means for magnetizing a volume an object in a first direction, the volume 
comprising a plurality of voxels; 

means for causing nuclear magnetic resonance magnetization to rotate into a 
plane transverse to the first direction; 

means for detecting the nuclear magnetic resonance field rotated into a transverse 
plane external to the object over a sample volume comprising a plurality of spatial 
locations; 

processor means for determining a set of components of spatial variation of the 
nuclear magnetic resonance field over a sample volume by association of the 
components by phase; 

processor means for determining a location of each voxel which is a source of the 
radio frequency field through each set of components representative of spatial 
variations of an intensity of the nuclear magnetic resonance field over the sample 
volume, and 
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processor means for displaying the voxel locations representative of anatomical 
images of the object based on selected physiological parameters. 

82. An apparatus for providing a two dimensional or three-dimensional nuclear 
5 magnetic resonance image of an object comprising: 

means for generating an image from a three-dimensional map of a transverse 
resonant radio frequency magnetic flux external to an object comprising a volume 
to be imaged, the volume comprising a plurality of voxels; and 

processor means for determining locations of voxels which is a source of a radio 
1 0 frequency field through spatial variations of an intensity of a nuclear magnetic 

resonance field of a set of components associated by phase. 

83. A method for providing a multidimensional image of inanimate and animate objects 
comprising the steps of: 

15 placing an object having a volume in a magnetic field, the volume comprising a 

plurality of voxels; 

applying a first source of radiation to the object; 

recording a secondary source of radiation emanating from the object at a plurality 
of spatial locations; 

20 forming components of the spatial variation of the secondary source of radiation 

external to the object due to voxels of the object; 

determining the location of voxels from the intensity variation of the components; 
generating and displaying the image from the locations and the intensity of the 
secondary radiation of each voxel. 

25 

84. The method according to claim 83, characterized in that the secondary recorded 
radiation is a three-dimensional transverse resonant radio frequency magnetic flux 
external to the object recorded as a function of time. 

30 85. The method according to claim 84, characterized in that the three-dimensional 

transverse resonant radio frequency magnetic flux external to the object is nuclear 
magnetic resonance radiation. 
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86. The method according to any one of claims 84 or 85 , further comprising recording 
the secondary source of radiation emanating from the object at a plurality of spatial 
locations using a plurality of detectors at each spatial location. 

5 87. The method according to any one of claims 83-86, characterized in that the recorded 
signals from each voxel contains components each having a unique phase angle 
relative to other components from different voxels. 

88. The method according to any one of claims 83-87, characterized in that a time 
10 dependent signal at each detector may be transformed into a series of components 

having intensity and phase data. 

89. The method according to any one of claims 83-88, characterized in that each set of 
components of the nuclear magnetic resonance signal over a sample volume due to a 

15 given voxel is determined from phase data and detector positions. 

90. The method according to any one of claims 83-89, characterized in that an intensity 
variation of a transverse radio frequency field over the sample volume is used to 
determine the coordinate location of each voxel. 

20 

91 . The method according to any one of claims 83-90, characterized in that the radio 
frequency field is a near field which is a dipole that serves as a basis element to form 
a unique reconstruction. 

25 92. The method according to any one of claims 83-91, characterized in that each volume 
element is reconstructed independently in parallel with all other volume elements 
such that a scan time is no greater than a nuclear free induction decay time. 

93. The method according to any one of claims 83-92, characterized in that said object is 
30 at least one selected from the group consisting of a human body, an animal, and 

tissue. 

94. A method for providing a multidimensional image of inanimate and animate objects 
comprising the steps of: 
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magnetizing a volume of an object to be imaged, the volume comprising a 
plurality of voxels; 

applying a radiation field to the object to be imaged when disposed in the 
magnetic field to produce a secondary radiation field emanating from the object; 
5 producing output signals from a detector in response to said secondary radiation 

at a plurality of spatial locations outside of the object as a function of time; 

determining a plurality of Fourier components each having the same frequency, 
an intensity and a phase angle; 

associating the Fourier components due to each voxel of the object to be imaged 
10 by phase to form a set corresponding to the spatial variation of intensity of the 

secondary radiation due to each voxel at the plurality of spatial locations outside of 
the object; 

converting each set of components into a voxel location; and 
producing an image based on an intensity of the secondary radiation from each 
15 voxel, the location of each voxel, and plotting the superposition of results for each 

voxel. 

95. A method for providing a two dimensional or three-dimensional nuclear magnetic 
resonance image of an object comprising the steps of: 

20 generating an image from a three-dimensional map of a transverse resonant radio 

frequency magnetic flux external to a volume of an object being scanned, the 
volume comprising a plurality of voxels, characterized in that a nuclear magnetic 
resonance signal of each voxel at any given detector gives rise to a Fourier 
component with a unique phase angle relative to the Fourier component of any other 

25 voxel of the object at a detector; 

a set of Fourier components that correspond to a nuclear magnetic resonance 
signal of a given voxel over detectors is determined; 

an intensity variation of the transverse radio frequency field over space is used to 
determine a coordinate location of each voxel; and 

30 a superposition of results for each voxel is plotted to provide a total image. 

96. A method of providing a multidimensional image of inanimate and animate objects 
comprising the steps of: 
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aligning magnetic moments of nuclei of a volume of an object to be imaged using 
a primary field, the volume comprising a plurality of voxels; 

further aligning the magnetic moments by a radio frequency pulse or series of 
pulses; 

5 recording free induction decay signals; 

Fourier transforming time dependent nuclear magnetic resonance signals to give 
intensity and phase of each Fourier component, characterized in that nuclear 
magnetic resonance signal of each voxel at any given detector gives rise to a Fourier 
component with a unique phase angle relative to the Fourier component of any 
10 other voxel at that detector; 

determining a matrix of Fourier components that correspond to the nuclear 
magnetic resonance signal of a given voxel over the detectors; 
determining matrices for a plurality of voxels; 

determining a coordinate location of each voxel from measurements of spatial 
15 variations of a transverse radio frequency field of each corresponding matrix; and 

processing and displaying the position of each voxel with a representation of the 
intensity of the nuclear magnetic resonance signal emitted by each voxel. 

97. A method according to any one of claims 83-96, characterized in that a strength and 
20 duration of rotating VL 1 radio frequency field that is resonant with the protons of the 

magnetized volume and is oriented perpendicularly to a direction of the magnetizing 
field is applied such that final precession angle of magnetization is 90 °, <£ Hi = 90° , 

such that a radio frequency dipole is transverse to a primary magnetizing field and 
perpendicular to a radio frequency magnetic field detector. 

25 

98. The method according to any of claim 83-97, characterized in that nuclear magnetic 
resonance pulse sequences which provide signals for a T x or T 2 image are applied. 



30 



99. The method according to claim 98, characterized in that the nuclear magnetic 

resonance pulse sequences which provide the signals for a or T 2 image comprises 
applying a 90 ° pulse followed by a series of 180 ° pulses. 
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100. The method according to claim 99, characterized in that the step of applying a 90 
° pulse followed by a series of 180 ° pulses comprises the Carr-Purcell-Meiboom- 
Gill sequence. 



5 101. The method according to any one of claims 96-100, characterized in that the step 
of determining the matrix of Fourier components that correspond to the nuclear 
magnetic resonance signal of a given voxel over the detectors is achieved by using a 
first component having a phase angle and calculating the phase angle as a function of 
spatial position of a first detector relative to any other detector and identifying a 
10 component at each detector having the calculated phase angle. 

102. The method according to any one of claims 96-100, characterized in that each 
matrix of components associated by phase comprises the intensity variation over the 
sample volume of the radio frequency field of bulk magnetization M of each voxel. 

15 

103. The method according to claim 96, further comprising reconstructing a 
multidimensional image representation according to spatial variations of the detected 
radio frequency magnetic field using at least one of a matrix inversion algorithm, a 
reiterative algorithm or Fourier transform algorithm. 

20 

1 04. The method according to claim 1 03, characterized in that a spatial Fourier 
transform algorithm is used. 

105. The method according to claim 104, characterized in that the spatial Fourier 
25 transform algorithm is performed on each set of components over the detector array 

to map each bulk magnetization M corresponding to a voxel to a spatial location of 
the imaged object. 



106. The method according to claims 104 or 105, characterized in that the spatial 
30 Fourier transform comprises the steps of : 

recording a nuclear magnetic resonance signal at discrete points in the sampling 
space at position x, y, z as a function of time; 

designating each point, x, y, z, as an element in matrix A ; 
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Fourier transforming time dependent signals to give an intensity and phase of 
each component wherein the nuclear magnetic resonance signal of each voxel at any 
given detector gives rise to a Fourier component with a unique phase angle relative 
to the Fourier component of any other voxel of the object at that detector; 
5 determining a matrix of Fourier components that correspond to the nuclear 

magnetic resonance signal of a given voxel over the detectors; 

determining matrices A n for all of the voxels; 

determining a coordinate location of each voxel using the measurements of the 
spatial variations of the transverse radio frequency field of a given matrix; and 
10 plotting a superposition of results for each voxel to provide a total image. 

107. The method according to claim 106, characterized in that the step of determining 
the coordinate location of each voxel using the measurements of the spatial variations 
of the transverse radio frequency field of a given matrix comprises the steps of 

a. ) discrete Fourier transforming each matrix A n to obtain a matrix B n ; 

b. ) multiplying each element of each matrix B n by a corresponding inverse value 
or reciprocal value of the Fourier transform of a system function which gives the 
spatial variation of the radio frequency field emitted by each voxel as a function of 
distance from the voxel to each detector, evaluated at the same frequency as the 
element of the matrix A n to form matrix C n ; 

c. ) generating matrix D n by taking the discrete inverse Fourier transform of 
matrix C n ; 

d. ) multiplying each element of each matrix D n by the distance squared along the 
z-axis to which the element corresponds to generate the position of the bulk 
magnetization M of voxel n which corrects the limitation of the sample volume to 
z>0; 

e. ) repeating steps a-d for all matrices A n ; and 

f. ) plotting superposition of results for each voxel to provide a total image. 

30 108. The method according to claims 106 or 107, characterized in that at each 

sampled external spatial location, x, y, and z in Cartesian coordinates, a magnitude 
B of the transverse radio frequency magnetic flux emitted by the corresponding z- 
axis aligned magnetic moment m z of each voxel is given by 



15 



20 



25 
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. _ 2z 2 -x 2 -y 2 
(x +y +z ) 

where the magnetic moment corresponds to the voxel's bulk magnetization M . 



10 



109. The method according to any one of claims 106-108, further comprising the step 
of determining a matrix of Fourier components that correspond to the nuclear 
magnetic resonance signal of a given voxel over the detectors by using a first 
component having a phase angle and calculating the phase angle as a function of 
spatial position of a first detector relative to any other detector and identifying a 
component at each detector having the calculated phase angle. 

110. The method according to claim 109, characterized in that each matrix of 
components associated by phase comprises the intensity variation over the sample 
volume of the radio frequency field of the bulk magnetization M of each voxel. 



15 111. The method according to claims 106 or 107, characterized in that a system 
function as a function of the detector location relative to each voxel in Cartesian 
coordinates, x, y, z, is by 
2z 2 ~x 2 ~y 2 



B = m 



r 2 t 2 , 2%5/2 

(x + y +z ) 



20 112. The method according to any one of claims 106-1 1 1, further comprising the step 
of correcting a point spread of each reconstructed voxel by assigning one voxel 
above a certain threshold with the bulk magnetization M and the other voxels are 
assigned a zero value. 

25 113. The method according to any one of claims 106-1 12, further comprising the step 
of correcting a finite size of each detector by multiplying each numeric value of 
matrix B n by the corresponding matrix element of the inverse or reciprocal of the 
Fourier transform of a waveform representative of a detector geometry. 



30 



114. The method according to claim 106-1 13, characterized in that each matrix 
element of the inverse or reciprocal of the Fourier transform of a waveform 
representative of the detector geometry is given by: 
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1 

smc(k x ,k y ,k z ) ^ 

115. The method according to any one of claims 106, 107, or 1 1 1-1 14, characterized 
in that each volume element is reconstructed independently in parallel with all other 

5 volume elements such that a scan time is no greater than a nuclear free induction 
decay time. 

116. The method according to any one of claims 106, 107, or 1 1 1-1 15 wherein the 
Fourier transform of the system function is given by: 



10 




117. The method according to any one of claims 106-1 16, characterized in that each 
matrix element for the step of multiplying each element of each matrix D n by the 
distance squared along the z-axis to which the element corresponds to generate the 
15 position of the bulk magnetization M of voxel n which corrects the limitation of the 
sample volume to z > 0 is given by: 

20 118. The method according to claim 103, characterized in that the step of 

reconstructing the multidimensional image uses a matrix inversion reconstruction 
algorithm. 

119. The method according to claim 118, characterized in that the matrix inversion 
25 reconstruction algorithm comprises the steps of : 

using the geometric system function to determine the spatial intensity variation of 
the transverse radio frequency field over a detector array to form a corresponding 
matrix; 

inverting the corresponding matrix; and 
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multiplying the signal over the detector array by the inverted matrix to give voxel 
sources. 



120. The method according to claim 119, further comprising the steps of 
5 a.) determining a coefficient for each voxel mathematically or by calibration 

which when multiplied by the bulk magnetization M of each voxel is that voxel's 
contribution to the signal at a given detector with the corresponding unique phase at 
each detector; 

b.) repeating step a for every detector; 
10 c.) using the coefficients to determine a matrix which when multiplied by a 

column vector of the bulk magnetization M values of the voxels, gives the signals at 
the detectors; 

d. ) inverting the matrix and storing the result in memory; 

e. ) recording signals as a function of time over the detector array; 

15 f.) Fourier transforming the signals as a function of time to give the Fourier 

components each having an amplitude and a unique phase; and 

g,) multiplying the components by the inverse matrix to generate a bulk 
magnetization M map which is displayed. 

20 121. The method according to claim 120, characterized in that the coefficient for each 
voxel determined mathematically is given by 

where V[x,y,z] is a detector signal at Cartesian coordinate x, y, z, C Q is a constant 
characteristic of each detector, l X9 / 2 , and l 3 are the total lengths of the detector array 
25 in the x, y, and z directions, & 15 k 2 , and k 3 are spacing of elements of the detector 

array in the x, y, and z directions, and n 19 n 2 , and n 3 are integer indices of the 
elements of the detector array in the x, y, and z directions. 
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122. 



The method according to claim 118, further comprising the steps of 
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determining a matrix of Fourier components that correspond to the nuclear 
magnetic resonance signal of a given voxel over the detectors; 

determining corresponding matrices for a plurality of voxels; 

determining a coordinate location of each voxel from measurements of spatial 
5 variations of the transverse radio frequency field of a given matrix wherein each 

matrix of components associated by phase comprises the intensity variation over the 
sample volume of the radio frequency field of the bulk magnetization M of each 
voxel; and 

multiplying each matrix of components by the inverse matrix of system function 
10 to generate the bulk magnetization M map. 



123. The method according to claim 122, further comprising the steps of 

a.) correcting a point spread of a reconstructed voxel by assigning one voxel 
above a certain threshold with the bulk magnetization M ; 
15 b.) assigning other voxels with a zero value; 

c. ) repeating steps a-b for other voxels; and 

d. ) plotting and displaying the superposition of the separate maps for each voxel 
which comprises the total bulk magnetization map as a nuclear magnetic resonance 
image. 

20 

124. The method according to claim 103, characterized in that said multidimensional 
image is reconstructed using a reiterative reconstruction algorithm. 



125. The method according to claim 124, characterized in that the reiterative 
25 reconstruction algorithm comprises the steps of 

a. ) using the geometric system function to determine a system of linear equations 
which gives the intensity, spatial variation, and phase of radio frequency field over a 
sample volume wherein the system of linear equations gives the signal and phase 
from each voxel at each detector based on the bulk magnetization M value of each 

30 voxel and the position of the voxel relative to the detector; 

b. ) Fourier transforming the signal as a function of time to give Fourier 
components each having an amplitude and a unique phase; 

c. ) determining weighting coefficients based on the linear equations; 

d. ) determining coefficients for each voxel mathematically or by calibration; 
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e. ) estimating bulk magnetization M for each voxel; 

f. ) calculating signals at each phase at each detector; 

g. ) comparing a result of multiplying the bulk magnetization M value of each 
voxel times its weighting coefficient and its calculated phase at a given detector to a 

5 measured signal and phase; 

h. ) making a correction to M of each voxel to give a second, or recomputed, 
estimate for M of each voxel; 

L) computing the signal value from this second estimate and making corrections 
according to steps f-h; 
10 j.) repeating steps f-i until the correction for each reiteration approaches a 

predefined limit which serves to indicate that the reconstruction is within reasonable 
limits of error; 

k.) repeating steps e-j of the procedure for other detectors; 
and 

1 5 1.) plotting and displaying a bulk magnetization map. 

126. The method according to claim 125, wherein the coefficient for each voxel 
determined mathematically is given by 

+y* +U2* -k/2* Efz+^f -\y-n-kf -fc-TLkf 1 

iwi-iz z z I f y J / , ;V 

«3=o ^=-^2* n^-ijiktix-nfaj +\y-n 2 k 2 J + [z + ^fcj j I 
20 where V[x,y,z] is the detector signal at Cartesian coordinate x, y, z, C Q is a constant 
characteristic of each detector, l v / 2 , and l 3 are the total lengths of the detector array 
in the x, y, and z directions, k u k 2 , and k 3 are the spacing of the elements of the 
detector array in the x, y, and z directions, and n 15 n 2 , and n 3 are the integer indices of 
the elements of the detector array in the x, y, and z directions. 

25 

127. A method of providing a multidimensional image of inanimate and animate 
objects comprising the steps of : 

a. ) recording secondary radiation signals corresponding to a transverse radio 
frequency magnetic field in response to a first radiation applied to a magnetized 

30 object over a sample volume comprising a plurality of voxels; 

b. ) forming a matrix of signals at a plurality of spatial positions each 
corresponding to a detector element of a detector array; 
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c. ) transforming the matrix numeric values according to a Fourier transform to 
give an intensity and phase of each component wherein a nuclear magnetic 
resonance signal of each voxel at any given detector gives rise to a Fourier 
component with a unique phase angle relative to the Fourier component of any other 

5 voxel of the object at that detector; 

d. ) determining a matrix of Fourier components that correspond to the nuclear 
magnetic resonance signal of a given voxel over the detectors; 

e. ) multiplying each numeric value of the transform array by a value which is an 
inverse or reciprocal of the Fourier transform of a waveform corresponding to an 

1 0 operation performed in step b; 

f. ) multiplying each numeric value of the discrete spectrum matrix by a value 
corresponding to an inverse Fourier transform of a function representative of the 
steps of exciting and detecting and providing a system corrected matrix; 

g. ) inverse-transforming the system corrected matrix according to a 

1 5 multidimensional discrete inverse Fourier transform providing a voxel location; 

h. ) correcting each element in the transformed matrix for any asymmetries in the 
dimensions of the sample volume; 

i. ) repeating steps c-h for other matrices, each corresponding to a distinct voxel; 
and 

20 j.) superimposing the results from steps a-h, plotting and displaying the image. 



128. The method according to claim 127, characterized in that at each sampled 

external spatial location, x, y, and z in Cartesian coordinates, the magnitude B of the 
transverse radio frequency magnetic flux emitted by the corresponding z-axis aligned 
25 magnetic moment m z of each voxel is given by 

2z 2 -x 2 ~y 2 



B = in 



f 2 , 2 , 2^5/2 

(x + y +z ) 



where the magnetic moment corresponds to the voxel's bulk magnetization M . 



30 129. The method according to claims 127 or 128, characterized in that the step of 
determining the matrix of Fourier components that correspond to the nuclear 
magnetic resonance signal of a given voxel over the detectors is achieved by using a 
first component having a phase angle and calculating the phase angle as a function of 
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spatial position of a first detector relative to any other detector and identifying the 
component at each detector having the calculated phase angle. 

130. The method according to claim 129, characterized in that each matrix of 

5 components associated by phase comprises an intensity variation over the sample 
volume of the radio frequency field of the bulk magnetization M of each voxel. 

131. The method according to any one of claims 127-130, characterized in that a 
system function as a function of the detector location relative to each voxel in 

10 Cartesian coordinates, x, y, z 5 is given by 

2z 2 -x 2 ~y 2 



B = ?n 



/ 2 , 2 , 2x5/2 



132. The method according to any one of claims 127-13 1, further comprising the step 
of correcting a point spread of each reconstructed voxel by assigning a voxel above a 
1 5 certain threshold with the bulk magnetization M and other voxels are assigned a zero 
value. 



133. The method according to any one of claims 127-132, further comprising step e of 
correcting a finite size of each detector by multiplying each numeric value of the 

20 matrix formed in step d by the corresponding matrix element of the inverse of the 
Fourier transform of a waveform representative of a detector geometry. 

134. The method according to claim 133, characterized in that each matrix element of 

an inverse of the Fourier transform of a waveform representative of the detector 

25 geometry is given by formula 

1 

smc(k x ,k y ,k 2 ) 



30 



135. The method according to any one of claims 127-135, characterized in that each 
volume element is reconstructed independently in parallel with all other volume 
elements such that a scan time is no greater than a nuclear free induction decay time. 
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136. The method according to any one of claims 127-135, characterized in that the 
step f of multiplying each numeric value of the discrete spectrum matrix by a value 
corresponding to the reciprocal of the Fourier transform of a function representative 
of the steps of exciting and detecting and providing a system corrected matrix 
5 wherein the discrete values of the Fourier transform of the system function matrix are 
given by: 




137. The method according to any one of claims 127-136, characterized in that the 
10 step h of correcting each element in the transformed matrix for any asymmetries in 
dimensions of the sample volume wherein the discrete values of are given by: 




138. The method according to claim 103, further comprising the steps of using the 
15 matrix inversion algorithm or the reiterative algorithm in combination with the 

Fourier Transform Algorithm. 

139. The method according to claim 103, further comprising the steps of using the 
matrix inversion algorithm and the reiterative algorithm in combination with the 

20 Fourier Transform Algorithm. 

140. The method according to claim 103, further comprising using the Fourier 
Transform Algorithm as input for the first reiteration of a reiterative algorithm. 

25 141. The method according to claim 140, further comprising the steps of: 

using the Fourier Transform Algorithm to calculate a series of first bulk 
magnetization images for each matrix of components associated by phase; 

detecting at least one peaks or edge of each image by using a cut-off level for the 
bulk magnetization assigned to a voxel or the value of the derivative of the bulk 
30 magnetization as a function of voxel position, respectively; 
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setting the bulk magnetization of the voxel or voxels corresponding to the peak or 
edge equal to the peak bulk magnetization and assigning the remaining voxels a zero 
value; 

calculating a signal at each detector using a system function; 
5 comparing each calculated signal to the actual signal and determining a 

difference or error; 

projecting the error over the voxels using the Fourier Transform Algorithm, 
Matrix Inversion Algorithm, or Reiterative Algorithm; 

adding back projected errors to a preceding corrected series of bulk 
10 magnetization maps; and 

performing the steps reiteratively until it stabilizes. 

142. The method according to claim 96, characterized in that a T x image is produced 
by a applying at least one pulse sequence that inverts the magnetization and records 

15 the relaxation. 

143. The method according to claim 142, further comprising producing a T x image by 
the steps of: 

switching on radio frequency receivers to follow a decay following a nuclear 
20 excitation comprising a H t pulse; 

examining M z (t), a time dependent bulk magnetization in the direction of a 
primary field H 0 , at a time t 0 after an inverting pulse by applying another H t pulse 
equivalent to a rotation by 90° after waiting the time t Q following the initial 
inversion, 

25 characterized in that the 90° pulse puts the z magnetization M z (/ 0 ) into the 

transverse plane for observation, and 

changing the waiting time t Q to allow for observation of M z (t) at different times 
during relaxation. 

30 144. The method according to claim 96, characterized in that a T 2 image is produced 
by a applying at least one pulse sequence that flips a magnetization vector into a 
transverse plane and records a transverse relaxation by producing at least one spin- 
echo. 
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145. The method according to claim 144, further comprising producing a T 2 image by 
the steps of 

flipping the magnetization initially aligned along the z R -axis into the transverse 
5 plane by applying a Hj pulse; 

recovering the T 2 relaxation from a T* FID by applying a 180° Hj pulse at time 
t Q along the y R -axis where the T* relaxation occurs in time t 0 ; 

recording the T 2 relaxation after an additional time t 09 characterized in that the 
total time elapsed after the 90° H x pulse is 2t Q = t E which is the spin-echo time and 
10 the dephasing due to the static inhomogeneity of H 0 is canceled out, whereby the 

relaxation is due to processes that create T 2 relaxation; 

applying a plurality of 180° pulses to form a plurality of spin echoes; and 

determining the T 2 relaxation from the envelop of the maximum amplitude of the 
spin echoes. 

15 

146. The method according to claim 145, further comprising a pulse sequence to give 
the data for a T 2 image known as the Carr-Purcell-Meiboom-Gill sequence 
comprises applying a 90° pulse along the x R -axis followed by a series of 180° pulses 
along the y R -axis at times t 0 + 2nt 0 where n is an integer including zero. 

20 

147. A method of providing a multidimensional nuclear magnetic resonance image of 
an animate or inanimate object comprising the steps of: 

placing a volume of an object having nuclei in a static magnetic field H 0 , the 
volume comprising a plurality of voxels, and then subjected them to an additional 
25 rotating radio frequency field , which is synchronous with their precession, such 

that the bulk magnetization M of each voxel of the object precesses about H x and 
rotates away from the primary field H 0 by an angle <f> Hi ; 

measuring of the intensity of the radio frequency signal over time and space 
following T x and/or T 2 encoding pulses wherein the magnitude of M is a maximum 
30 initially and decays with time by emission of the same multipolarity radiation that it 

absorbed and by transfer of energy to the surrounding lattice and the intensity of the 
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radiation is a function of M and the coordinate position relative to the radio 
frequency emitting voxel; 

Fourier transforming the signal as a function of time at a given detector position 
to give the components each having an amplitude and a unique phase; 
5 determining each set of components of the nuclear magnetic resonance signal 

over the sample volume due to a given voxel from the phase data and the detector 
positions; 

determining the location of each voxel through the spatial variations of the 
intensity of the transverse nuclear magnetic resonance field of the set of components 
1 0 associated by phase; and 

superimposing the results of the determination of the location of M of each voxel 
and the plotting and displaying the resulting image. 



148. The method according to claim 96, characterized in that the plurality of Fourier 
15 components each having the same frequency, an intensity and a phase angle each 
have a unique phase. 



149. The method according to claim 148, characterized in that the unique phase of 
each component is provided by a combination of 
20 the angle 6 suspended between each spatial location at which the secondary field 

is detected and a radial vector from each voxel to each spatial location; and 

the angle 0 due to a separation distance r between each voxel and each spatial 
location given by a wavenumber of a detected field k times r . 



25 150. The method according to claim 87, characterized in that the radio frequency field 
is sampled synchronously so that the phase at any given detector may be related to 
that at any other detector. 

151. The method according to claim 87, characterized in that the radio frequency field 
30 is sampled at known times so that the phase at any given detector may be related to 
that at any other detector. 
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152. The method according to claim 96, characterized in that the processor which 
determines a plurality of Fourier components each having the same frequency, an 
intensity and a phase angle uses a first component having a phase angle and 
calculates the phase angle as a function of spatial position of the first spatial position 

5 relative to any other spatial location and identifies the component at each spatial 
location having the calculated phase angle. 

153. The method according to claim 152, characterized in that the voxels that are on 
the same axis with a given detector or a plurality of detectors that align vertically 

10 with one or more voxels is determined. 

154. The method according to claim 153, characterized in that the phase angle is a 
function of the distance r of the detector from the voxel. 

15 155. The method according to claim 154, characterized in that the components may 
be identified by the presence of at least one component of a phase given by 4- kAr 
where ^ is the phase angle of the first component, Ar is the distance of the second 
detector relative to the first of the detector, and k is the wavenumber of the nuclear 
magnetic resonance radiation emanating from the object wherein this detector(s) is 

20 defined as having a zero angle 9 . 

156. The method according to claim 155, characterized in that a phase angle <j> 2 of a 
second component at a second detector aligned on the detection axis corresponding 
to a first component at a first detector having a phase angle ^ is given by 
25 $ 2 = 0! + Mr 

where Ar is the distance of the second detector relative to the first of the detector 
and k is the wavenumber of the nuclear magnetic resonance radiation emanating 
from the object. 

30 157. The method according to claim 156, characterized in that, from the components 
identified as coaxial with a detector, the phases of the components from the voxel are 
determined at the other detectors using the relative position of other detectors of the 
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sample volume relative to the first detector and <j> x of the component at the first 
detector. 



158, The method according to claim 157, characterized in that the phase 0 due to a 
5 separation distance between a voxel and a detector of r is given by the wavenumber 

of the radio frequency field k times r 9 
(j) =kr 

where the wavenumber k is given by 

2k 
~ % 

10 the phase angel (j) is given by 

2tc 



where Ax, Ay, and Az is the separation distance between a voxel and a detector in 
the x, y 9 and z directions, respectively, and 

in the case that a two dimensional slice is imaged in the transverse plane, the 
15 phase angel <j> is given by 

^=*r = ^V(A*) 2 + (Ay) 2 



159. The method according to claim 158, characterized in that all voxels which have 
the same radial distance from a given detector will have the same phase (j) . 

20 

160. The method according to claim 159, characterized in that at each case that two 
detectors have the same phase (p, since the radial distances r of two voxels to a 
given detector are equivalent, a unique angle 6 suspended between the direction of 
the detector and the vector along r exists. 

25 

161. The method according to claim 160, characterized in that the radio frequency 
dipoles of all voxels are time synchronous. 



162. The method according to claim 96, characterized in that the nuclear magnetic 
30 resonance pulsing and detection may be synchronized. 
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163. The method according to claim 96, characterized in that at least one standard 
dipole may be used as reference phase to set the phases of the radio frequency 
dipoles of the voxels. 

5 164. The method according to claim 163, characterized in that all voxels that are not 
collinear with the axis of a detector posses an angle 0 suspended between the 
direction of the detector and the radial vector from the dipole to the detector. 

165. The method according to claim 164, characterized in that since the radio 

10 frequency dipoles are time synchronous, at each point in time, the angle corresponds 
to a phase angle 0 of the radio frequency signal from each voxel at each detector. 

166. The method according to claim 165, characterized in that in the case that a two 
dimensional slice is imaged in the transverse plane, the phase angle 0 is given by 



167. The method according to claim 166, characterized in that the total unique phase 
angle (j) T of the radio frequency signal from each voxel at each detector is given as 
the sum of (j) and 9 given by 



15 
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<j> T =<f>+0 
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